求平方根算法 Heron’s algorithm

Posted 2022-08-21 data_algorithms

求平方根问题

• 概述：本文介绍一个古老但是高效的求平方根的算法及其python实现，分析它为什么可以快速求解，并说明它为何就是牛顿迭代法的特例。

• 问题：求一个正实数的平方根。

  给定正实数 \(m\)，如何求其平方根\(\sqrtm\)？
  你可能记住了一些完全平方数的平方根，比如\(4, 9, 16, 144, 256\)等。那其它数字（比如\(105.6\)）的平方根怎么求呢？

  实际上，正实数的平方根有很多算法可以求。这里介绍一个最早可能是巴比伦人发现的算法，叫做Heron’s algorithm。

算法描述

假设要求的是\(x = \sqrtm\)，Heron’s 算法是一个迭代方法。
在迭代的第\(k=0\)步，算法随机找一个正数作为\(x\)的初始值，记为\(x_0\)，例如\(x_0 = 1\)。
在第\(k \in [1, 2, ...]\)步，计算 \(x_k = \fracx_k-12 + \fracm2 x_k-1\)。
每迭代一步，算法计算\(x_k\)与\(x_k-1\)的变化，\(\Delta_k = |x_k - x_k-1|\)，若 \(\Delta_k < \epsilon\)，则算法结束，输出\(x_k\)。
这里\(\epsilon\)是计算精度要求，可以取一个小正数，如\(1e-14\)。算法描述如下：

1. 初始化 \(x = 1\);
2. 计算 \(\Delta_k = |x_k - x_k-1|\);
3. 若\(\Delta_k < \epsilon\)，返回\(x\)，算法结束;
4. \(x \leftarrow \fracx2 + \fracm2 x\);
5. 返回第2步;

在 \(m\)比较大时，可能会出现迭代次数增多，数值不稳定的情况。
我们可以通过将\(m\)缩放到一个较小的区间，求缩放后的平方根，然后再反缩放得到\(m\)的平方根。
考虑任意一个正实数 \(m\)，我们总是可以将\(m\)看作两个正实数的乘积 \(m = n * 4^u\)，其中\(n \in [0.5, 2]\)。
此时有 \(\sqrtm = \sqrtn * \sqrt2^2u = 2^u * \sqrtn\)。
因此只需要计算 \(\sqrtn\)，计算结果乘以\(2^u\)，就可以得到 \(\sqrtm = 2^u *\sqrtn\)。
由于\(n \in [0.5, 2]\)，求根过程更容易，误差也可以得到较好的控制。
这里的\(u\)是将\(m\)连续\(u\)次除以4，或者将\(m\)连续\(u\)次乘以4来缩放到指定的\([0.5, 2]\)区间：

1. 初始化 \(u=0\);
2. 若 \(0.5 \leq m \leq 2\),返回 \(u\)，算法结束；
3. 若 \(m > 2\): \(m \leftarrow m / 4\), \(u \leftarrow u+1\);
4. 若 \(m < 1\): \(m \leftarrow m * 4\), \(u \leftarrow u-1\);
5. 返回第2步;

代码实现

"""
Heron's 求平方根算法
@data_algorithms
"""


def heron(m, epsilon=1e-15):
    """
    Heron's 求根算法
    @m: 带求根的正实数
    @epsilon: 迭代结束条件
    @return: m的平方根
    """
    if m <= 0:
        raise ValueError("Non-negative real numbers only.")
    m, u = scale(m) #m缩放到[0.5, 2]区间
    x = 1
    delta = abs(x)
    while delta >= epsilon:
        newx = 0.5 * (x + m / x)
        delta = abs(newx - x)
        x = newx
    return x * (2 ** u)


def scale(m, base=4):
    """
    @m: 正实数m
    @base: m = base ^ u + m~
    @return: m~, u
    """
    u = 0
    while m > 2:
        m = m / 4
        u += 1
    while m < 0.5:
        m = m * 4
        u -= 1
    return m, u


def heron_test():
    from math import sqrt
    from random import random
    epsilon = 1e-15
    for _ in range(100000):
        m = random()* 1e10
        error = abs(heron(m) - sqrt(m))
        assert error < 1e-10

if __name__ == "__main__":
    from math import sqrt
    for x in [105.6, 0.1, 0.5, 1.5, 2, 
              4, 10, 123, 256, 1234]:
        print(x, heron(x), sqrt(x))
    heron_test()

算法分析

为什么Heron’s 算法能够快速找到平方根呢？

我们可以通过观察每一步迭代的误差来进行分析。

假设要求解的真值为 \(x\)，即\(x=\sqrtm \Rightarrow m = x^2 = m\)。
在第\(k\)步，算法的误差是\(e_k = (x_k - x)\)。
在第\(k+1\)步，\(x_k+1 = \fracx_k2 + \fracm2 x_k\)， 其误差是 $e_k+1 = (x_k+1 - x) =(\fracx_k2 + \fracm2 x_k - x) = (\fracx_k - x2 + \fracm - x_k x2x_k) =(\frace_k2 - \fracx(x_k - x)2x_k) = (\frace_k2 - \fracxe_k2x_k) = \frace_k^22x_k $。
所以 \(e_k+1 = \frace_k^22x_k\)，即后一次迭代的误差与前一次的平方成正比。
只要某一步的误差在\((0, 1)\)区间内，则误差会快速的缩小，所以算法可以快速地找到平方根。

与牛顿法的关系

这个算法其实就是在用牛顿法求一个函数的根，所以是牛顿法的一个特列。

为什么呢？牛顿法求根的迭代公式是：\(x_k = x_k-1 - f(x_k-1) / f'(x_k-1)\)。
这里令 \(f(x) = m - x^2\)，则有 \(f'(x) = -2x\)，所以 \(f(x) / f'(x) = -(m - x^2)/(2x)\)，牛顿法迭代公式就变为 \(x_k = x_k-1 + (m - x_k-1^2)/(2x_k-1) = x_k-1 + m / (2x_k-1) - x_k-1 / 2 = x_k-1 / 2 + m / (2x_k-1)\)，这也就是Heron‘s的迭代公式了。