长除法计算平方根的方法总结与代码实现(C++, Python)

---

tags: DSA Math C++ Python

写在前面

之前总结了计算平方根的方法, 但是并没有给出手算方法的解释, 这次专门写一下手算方法.

据说这个方法是中国的数学家创造的, 我也没深入考证过, 总之就是非常经典了, 因为这个长除法算法(英文:Long Division Algorithm)可以计算任意精度的平方根, 也就是可以算小数点后的任意位, 下面来看看具体的方法与原理.

原理解释

代数

其实原理是基于这样一个式子:
$$x^2=(10a+b{)}^2⟺x^2-100a^2=(20a+b)b.\text{\hspace{1em}(*)}$$
就是说对于一个两位数$x$, 其平方(设其有4位)有这样的一种表示, 那么如果要计算某一个数$y=x^2$的平方根, 只需要通过长除法, 根据数的前面两位和后面两位迭代计算即可.

当然这样直接说显得有点不够直观, 我们举个例子, 对于$y=6561$, 有
$$\begin{array}{rl}& 8,1\\ 8& \sqrt{65,61}\\ & 64\\ 16\underline{1}& \overline{1,61}\\ & \underline{1,61}\\ & 0\end{array}$$

• 首先找到前两位(不妨设为$x_1$, $x_1=65$)的小于等于该两位数的最大整数$k$($k\in [0,9]$), 该数满足$k^2⩽x_1$, 那么显然有$k^2=8^2=64⩽65$, 这步之后, 其实就找到了商$a$, $a=k=8$.
• 然后计算余数, 即上面公式$(\ast )$中的$x^2-100a^2$, 这个值等于$161$(长除法中表现为借位),
• 最后去找数字$b$使得$(20a+b)b=(160+b)b⩽161$的最大的$b$, 其中$a=8$, 也就是上面找到的商. 显然$b=1$, 可以恰好整除余数.

上面的例子中给出的是恰好整除(平方根为正整数)的情况, 那么对于不能整除的情况呢?

对于平方根为无理数的情况, 上面式子$(\ast )$仍然成立, 只不过对应为余数始终不为$0$, 这样一直做, 就能得到平方根了.

几何

在YouTube看到一个很不错的对于长除法的原理解释, 通过分割正方形的方法来给出直观的几何解释, 大家可以看看, 我传了B站. 長除法開方的原理（粵語中文字幕）_哔哩哔哩_bilibili;

C++实现

代码方面一开始我不太熟悉, 参考了¹, 后来自己想出来了一种基于二分法的方法, 在寻找商数的时候比¹的代码简洁一些, 不用遍历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, 得到的值简直完美, 这里对照了oeis²给出的结果:

2.236067977499789696409173668731276235440618359611525724270897245410520

以及Python自带的求平方根函数的结果:

from math import sqrt
print('%.100f'%sqrt(5))

2.2360679774997898050514777423813939094543457031250000000000000000000000000000000000000000000000000000

发现Python自带的数值计算在20位之后就会出现误差了, 但是使用长除法就不会有误差.

ref

---

1. Long Division Method to find Square root with Examples - GeeksforGeeks; ↩︎ ↩︎

2. A002163 - OEIS; ↩︎