求平方根算法 Heron’s algorithm

求平方根问题

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

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

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

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

算法描述

假设要求的是\(x = \sqrt{m}\)，Heron’s 算法是一个迭代方法。
在迭代的第\(k=0\)步，算法随机找一个正数作为\(x\)的初始值，记为\(x_0\)，例如\(x_0 = 1\)。在第\(k \in [1, 2, ...]\)步，计算 \(x_{k} = \frac{x_{k-1}}{2} + \frac{m}{2 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 \frac{x}{2} + \frac{m}{2 x}\);
5. 返回第2步;

在 \(m\)比较大时，可能会出现迭代次数增多，数值不稳定的情况。我们可以通过将\(m\)缩放到一个较小的区间，求缩放后的平方根，然后再反缩放得到\(m\)的平方根。
考虑任意一个正实数 \(m\)，我们总是可以将\(m\)看作两个正实数的乘积 \(m = n * 4^u\)，其中\(n \in [0.5, 2]\)。
此时有 \(\sqrt{m} = \sqrt{n} * \sqrt{2^{2u}} = 2^u * \sqrt{n}\)。因此只需要计算 \(\sqrt{n}\)，计算结果乘以\(2^u\)，就可以得到 \(\sqrt{m} = 2^u *\sqrt{n}\)。
由于\(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=\sqrt{m} \Rightarrow m = x^2 = m\)。
在第\(k\)步，算法的误差是\(e_k = (x_k - x)\)。在第\(k+1\)步，\(x_{k+1} = \frac{x_{k}}{2} + \frac{m}{2 x_{k}}\)， 其误差是 $e_{k+1} = (x_{k+1} - x) =(\frac{x_{k}}{2} + \frac{m}{2 x_{k}} - x) = (\frac{x_k - x}{2} + \frac{m - x_k x}{2x_k}) =(\frac{e_k}{2} - \frac{x(x_k - x)}{2x_k}) = (\frac{e_k}{2} - \frac{xe_k}{2x_k}) = \frac{e_k^2}{2x_k} $。所以 \(e_{k+1} = \frac{e_k^2}{2x_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的迭代公式了。

原文：https://www.cnblogs.com/data_algorithms/p/heron_algorithm_square_root_python.html