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

 
Category: DSA

写在前面

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

据说这个方法是中国的数学家创造的, 我也没深入考证过, 总之就是非常经典了, 因为这个长除法算法(英文: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;

截屏2022-12-23 21.56.50

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位之后就会出现误差了, 但是使用长除法就不会有误差.

ref