写在前面
之前总结了计算平方根的方法, 但是并没有给出手算方法的解释, 这次专门写一下手算方法.
据说这个方法是中国的数学家创造的, 我也没深入考证过, 总之就是非常经典了, 因为这个长除法算法(英文:Long Division Algorithm)可以计算任意精度的平方根, 也就是可以算小数点后的任意位, 下面来看看具体的方法与原理.
原理解释
代数
其实原理是基于这样一个式子: \(x^2=(10a+b)^2\iff x^2-100a^2=(20a+b)b.\tag{*}\) 就是说对于一个两位数$x$, 其平方(设其有4位)有这样的一种表示, 那么如果要计算某一个数$y=x^2$的平方根, 只需要通过长除法, 根据数的前面两位和后面两位迭代计算即可.
当然这样直接说显得有点不够直观, 我们举个例子, 对于$y=6561$, 有 \(\begin{aligned} &\ \quad8,\ \ 1\\ 8&\sqrt{65,61}\\ & \quad64\\ 16\underline{1}&\quad\overline{\ \ 1,61}\\ &\quad\ \ \underline{1,61\,}\\ &\qquad\quad0 \end{aligned}\)
- 首先找到前两位(不妨设为$x_1$, $x_1=65$)的小于等于该两位数的最大整数$k$($k\in[0,9]$), 该数满足$k^2\leqslant x_1$, 那么显然有$k^2=8^2=64\leqslant65$, 这步之后, 其实就找到了商$a$, $a=k=8$.
- 然后计算余数, 即上面公式$(*)$中的$x^2-100a^2$, 这个值等于$161$(长除法中表现为借位),
- 最后去找数字$b$使得$(20a+b)b=(160+b)b\leqslant 161$的最大的$b$, 其中$a=8$, 也就是上面找到的商. 显然$b=1$, 可以恰好整除余数.
上面的例子中给出的是恰好整除(平方根为正整数)的情况, 那么对于不能整除的情况呢?
对于平方根为无理数的情况, 上面式子$(*)$仍然成立, 只不过对应为余数始终不为$0$, 这样一直做, 就能得到平方根了.
几何
在YouTube看到一个很不错的对于长除法的原理解释, 通过分割正方形的方法来给出直观的几何解释, 大家可以看看, 我传了B站. 長除法開方的原理（粵語中文字幕）_哔哩哔哩_bilibili;
C++实现
代码方面一开始我不太熟悉, 参考了1, 后来自己想出来了一种基于二分法的方法, 在寻找商数的时候比1的代码简洁一些, 不用遍历0~9, 效率也相对高一些.
#include <iostream>
#include <vector>
#include <iomanip>
#include <typeinfo>
using namespace std;
ostream& operator<<(ostream& os, const vector<int> v) {
for (int i : v) os << i << " ";
return os << endl;
}
int find_nice(int R, int b = 0) {
int l{}, r{9};
while (l <= r) {
int mid = l + (r - l) / 2;
if ((20 * b + mid) * mid > R)
r = mid - 1;
else
l = mid + 1;
}
return l - 1;
}
long long mySqrt(long long n) {
long long dividend{}, quotient{}, reminder{};
int i{};
vector<int> a(15, 0), quot{};
// Dividing the number into segments
while (n) {
a[i++] = n % 100;
n /= 100;
}
// i = 10;
// a[i - 1] = 5;
cout << a;
for (int j = i - 1; j >= 0; --j) {
dividend = reminder * 100 + a[j]; // update dividend
long long tmp = find_nice(dividend, quotient);
quot.emplace_back(tmp);
reminder = dividend - (20 * quotient + tmp) * tmp;
quotient = quotient * 10 + tmp;
// cout << quotient << typeid(tmp).name() << endl;
}
cout << quot;
return quotient;
}
void t1() {
cout << mySqrt(500) << endl; // 22
cout << mySqrt(839) << endl; // 28
cout << mySqrt(1009) << endl; // 31
}
void t2() {
cout << mySqrt(500000000000000L) << endl;
/* 2 2 3 6 0 6 7 9
22360679
*/}
int main(int argc, char const* argv[]) {
// t1();
t2();
return 0;
}
使用C++实现并不复杂, 但是却因为数值类型的位数要求, 导致结果总是会有误差的. 联想到Python强大的任意精度计算, 决定用Python来实现.
Python实现
用Python重写上面的代码, 可以得到任意精度的平方根值.
def find_nice(R, b=0):
l, r = 0, 9
while l <= r:
mid = l + (r - l) // 2
if (20 * b + mid) * mid > R:
r = mid - 1
else:
l = mid + 1
return l - 1
def mySqrt(n=0):
dividend = quotient = reminder = 0
a = [0] * 150
# quot = []
i = 150
a[i - 1] = 5
for j in range(i - 1, -1, -1):
dividend = reminder * 100 + a[j]
tmp = find_nice(dividend, quotient)
# quot.append(tmp)
reminder = dividend - (20 * quotient + tmp) * tmp
quotient = quotient * 10 + tmp
# print(quot)
return quotient
if __name__ == '__main__':
# print('%100.100f' % 5**.5)
print(mySqrt())
得到的结果如下:
223606797749978969640917366873127623544061835961152572427089724541052092563780489941441440837878227496950817615077378350425326772444707386358636012153
不得不说Python的任意精度数才是数值计算的不二之选, 利用Python重写上面的程序, 重新计算根号5, 得到的值简直完美, 这里对照了oeis2给出的结果:
2.236067977499789696409173668731276235440618359611525724270897245410520
以及Python自带的求平方根函数的结果:
from math import sqrt
print('%.100f'%sqrt(5))
2.2360679774997898050514777423813939094543457031250000000000000000000000000000000000000000000000000000
发现Python自带的数值计算在20位之后就会出现误差了, 但是使用长除法就不会有误差.