如何笔算开平方根

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

分块

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

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

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

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

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

初商和次商分别是什么

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

这个时候“？”的值应该是 2，那么此时“二十倍初商加次商”的值就是 $20 \times 1 + 2 = 22$ ，然后 $22 \times 2 = 44 < 67$ ，非常合理。

此时我注意到：

(20 \times 1 + 3) \times 3 = 69

而 $169 = 13^{2}$ ，于是我明白了 $20$ 这个常数大概跟十进制与平方和公式有关。

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

(10 a + b)^{2} = 100 a + 20 a b + b^{2} = 100 a + (20 a + b) b

证毕。

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

然后计算器算一下答案是 12.9228，算法正确，如果是拿着上一个结果 $2 \times 20 = 40$ ，那么即使是 $49 \times 9$ 也远远小于 $2300$ ，因此猜想正确。

编写程序

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

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$ 的模拟：

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

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

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

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

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)

稍微算一下复杂度，设这个数为 $m$ ，_20_p1 的这个过程是与 $m$ 的位数有关，即 $O (\log m)$ ，而每次 iter 也是算一位，假设算到整数位就停止，复杂度也是 $O (\log m)$ 。因此总的时间复杂度是 $O (\log^{2} m)$ 。

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

完整代码

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

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

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))