如何笔算开平方根

发表于 更新于 分类于 阅读次数： 本文字数： 3.9k 阅读时长 ≈ 4 分钟

小学的时候家父给我炫耀他上学的时候学过手算开平方根，说以后在朋友间炫耀可能有用。当时我可能学会了，但之后这个方法似乎在数学教材中绝版了，只记得“二十倍初商加次商”这个口诀。这篇文章的主要目的就是叙述，我如何回忆起来的过程。

分块

首先，我挑了一个 167 这个数进行测试：

稍微估计一下就可以知道答案在 12 到 13 之间，但竖式的上面有三个空格，而商却只有两位，所以得构思一下这两个数字应该写在什么位置。

换句话说，是应该用 16 去进行一次“二十倍初商加次商”这个操作，还是用 1 去操作。

显然，第一个数字 1 应该写在最高位上面，但这对所有的数都成立吗？于是我把 167 改成了 67：

答案在 8 到 9 之间，这个时候，如果对高位 6 进行操作，得到的结果就会是 2 开头，那显然就不对了。因此，我估计大概是以个位为基准，两位两位的进行操作。（因为一位数的平方可以“占据”两位数的空间，两位数的平方可以“占据”四位数的空间，因此两个两个一组是合理的）

初商和次商分别是什么

我们回到之前的数 167，首先对百位进行开根，我们有：

这个时候“？”的值应该是 2，那么此时“二十倍初商加次商”的值就是 ，然后 ，非常合理。

此时我注意到：

而 ，于是我明白了 这个常数大概跟十进制与平方和公式有关。

于是我们可以证明算法的正确性，设初商为，次商为，则：

证毕。

这便同时说明了“二十倍初商加次商”和两个两个一组的原理。此时次商很明显就是我们正在算的那一位，然后我们可以猜想初商应该是先前算出来的商的数值，于是我们接着算下去：

然后计算器算一下答案是 12.9228，算法正确，如果是拿着上一个结果，那么即使是 也远远小于 ，因此猜想正确。

编写程序

显然，手动计算验证费时费力，如果可以通过编写程序来加快用例验证，可以更加增强自己对算法的信心。于是我们首先编写“二十倍初商加次商”的函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def _20_p1(first, rem, next):
    """
    recursively tries largest x from 0 to 9 to get
    (20 * first + x) * x <= rem * 100 + next
    return x and remainer
    """
    x = 0

    for i in range(11):
        if (20 * first + i) * i > rem * 100 + next:
            x = i - 1
            break
        assert i != 10
    
    return (x, rem * 100 + next - (20 * first + x) * x)

这里的 first 指“初商”；rem 指当前遗留的余数，例如上一个例子的“0”，“23”，“59”；而 next 指下一位的数字，例如 “67”，“0”。

然后我们尝试手动使用 _20_p1 来模拟之前的算法逻辑，第一步得到 这个显然步骤我们暂时省略，直接来到得到第二位 的模拟：

1
2
3
prev, modulo = 1, 0
prev, modulo = prev * 10 + _20_p1(prev, modulo, 67)[0], _20_p1(prev, modulo, 67)[1]
print(_20_p1(prev, modulo, 0))

先前我们得到了“初商” prev 为 ，而 rem（即代码中的 modulo）为 ，于是我们使用 _20_p1(prev, modulo, 67) 得到了这一轮“二十倍初商加次商”的“次商”与余数，分别用后缀[0]与[1]表示。

这一轮的“初商”通过 合并现在的“次商”成为了下一轮的“初商”，而余数又成为新一轮的 modulo 变量被下一轮的 _20_p1(prev, modulo, 00) 加上 00，开始新一轮的迭代。

最后一个问题就是，如何获取next这个变量呢？答案很简单，把被除数看作一个字符串，两个两个分块，之后转回去即可。

于是，一个递推过程就这么形成了，边界条件（最开始那个开根）也很简单，直接将初商和 modulo 都设为 即可：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
res = ''
chunks = split_two_by_two(raw_number, int_digit % 2)

# now do the algorithm
prev, modulo = 0, 0
for i in range(iters):
    if i < len(chunks):
        next_chunk = int(chunks[i])
    else:
        next_chunk = 0

    x, rem = _20_p1(prev, modulo, next_chunk)
    prev, modulo = prev * 10 + x, rem
    res = str(prev)

稍微算一下复杂度，设这个数为，_20_p1 的这个过程是与 的位数有关，即 ，而每次 iter 也是算一位，假设算到整数位就停止，复杂度也是 。因此总的时间复杂度是 。

这个时间复杂度甚至不如用二分法算，时间复杂度为 ，所以“二十倍初商加次商”也只能在朋友面前炫耀一下了。

完整代码

通过了一些简单的 debug 和特判处理，我这里写出了具有一定健壮性的代码，通过了一系列测试用例。因此这里分享如下：

这里输入只能是自然数，还没有对小数提供支持（虽然也只是一大堆细节处理）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def _20_p1(first, rem, next):
    """
    recursively tries largest x from 0 to 9 to get
    (20 * first + x) * x <= rem * 100 + next
    return x and remainer
    """
    x = 0

    for i in range(11):
        if (20 * first + i) * i > rem * 100 + next:
            x = i - 1
            break
        assert i != 10
    
    return (x, rem * 100 + next - (20 * first + x) * x)

def getintdigit(x):
    return len(str(x))

def split_two_by_two(s, first_len):
    if first_len == 0:
        first_len += 2
    if first_len >= len(s):
        return [s]

    first_chunk = s[:first_len]
    remaining = s[first_len:]
    
    rest_chunks = [remaining[i:i+2] for i in range(0, len(remaining), 2)]
    
    return [first_chunk] + rest_chunks

def mysqrt(x, iters):
    """
    Returns the square root of x.
    iters means keep {iters} valid digits.
    """
    if not isinstance(iters, int) or (iters < 1):
        raise ValueError('iters must greater than 0')
    if not isinstance(x, int):
        raise ValueError('x must be an integer')
    if x < 0:
        raise ValueError('x must be non-negative')

    int_digit = getintdigit(x)
    res_int_digit = (int_digit - 1) // 2 + 1

    raw_number = str(x).replace(".", "")
    res = ''
    chunks = split_two_by_two(raw_number, int_digit % 2)

    # now do the algorithm
    prev, modulo = 0, 0
    for i in range(iters):
        if i < len(chunks):
            next_chunk = int(chunks[i])
        else:
            next_chunk = 0

        x, rem = _20_p1(prev, modulo, next_chunk)
        prev, modulo = prev * 10 + x, rem
        res = str(prev)

    if iters <= res_int_digit:
        return res + '0' * (res_int_digit - iters)
    return res[:res_int_digit] + '.' + res[res_int_digit:]

print(mysqrt(167, 10))