8.3 Heron算法

巴比伦算法

  求根号可以说是刚需了，最简单的就是计算$\sqrt{2}$，小时候做数学题我们是背下来了，就是等于$1.414$，但是三位小数是肯定精度不够的，需要再精确怎么办？哉4000年前，巴比伦人发明了一种求$\sqrt{2}$的算法：
$$\sqrt{2}=1+\frac{24}{60}+\frac{51}{60^2}+\frac{10}{60^3}$$
  那么这到底有多精确呢？我们用python试一试：

import math


def sqrt2():
    return 1 + 24 / 60 + 51 / (60 * 60) + 10 / (60 * 60 * 60)


if __name__ == '__main__':
    a = sqrt2()
    print(a, a * a)
    b = math.sqrt(2)
    print(b, b * b)

  执行结果：

1.414212962962963 1.9999983046124827
1.4142135623730951 2.0000000000000004

  精度还是不错的，但是远不如python自带的math库精确。但是在4000年前，巴比伦算法还是很实用的，只是这个时代，它的精度不够看了。

牛顿-拉普森-辛普森算法

  这个算法，名字好长啊，英文是Newton-Raphson-Simpson method，算法很简单，就是无限迭代。假设要开方的数字是$y$，开方后的数值是$x$，先将$y$赋值给x，然后不断迭代以下过程：
$$x\leftarrow \frac{x+y/x}{2}$$
  这个算法有多精确呢，我们用python迭代十次试试：

import math


def sqrt(y):
    x = y
    for i in range(10):
        x = (x + y / x) / 2
    return x


if __name__ == '__main__':
    a = sqrt(3)
    print(a, a * a)
    b = math.sqrt(3)
    print(b, b * b)

  执行结果：

1.7320508075688772 2.9999999999999996
1.7320508075688772 2.9999999999999996

  哇哦，不得了不得了，跟math库结果一模一样的。

Heron算法

  牛顿-拉普森-辛普森算法就是Heron算法的现代名称。但是这个算法怎么证明呢？我们可以将上述迭代十次的实验代码写得更正式一点：

import math


def sqrt(y):
    x = y
    while True:
        t = (x + y / x) / 2
        if t == x:
            break
        else:
            x = t
    return x


if __name__ == '__main__':
    y = 888
    a = sqrt(y)
    print(a, a * a)
    b = math.sqrt(y)
    print(b, b * b)

  执行结果：

29.79932885150268 888.0
29.79932885150268 888.0

  代码改成了迭代后x不再变化时，就跳出迭代，认为找到了最精确的值，受限于浮点数编码，这个迭代肯定会结束的，不会无限循环。在数学上看，这个迭代永远不会结束，在循环无限次后，才会使得x不再变化，也就是：
$$\begin{gathered} x=\frac{x+y/x}{2} \\ \Rightarrow x^2=\frac{x^2+y}{2} \\ \Rightarrow \frac{x^2}{2}=\frac{y}{2} \\ \Rightarrow x^2=y \end{gathered}$$
  所以x的初始值是不重要的，只要迭代次数多，最终都会得到$x=\sqrt{y}$，我们可以大胆尝试下，初始将$x$赋值为一个随机数，试一下：

import math
import random


def sqrt(y):
    x = random.randint(1,10000)
    while True:
        t = (x + y / x) / 2
        if t == x:
            break
        else:
            x = t
    return x


if __name__ == '__main__':
    y = 888
    a = sqrt(y)
    print(a, a * a)
    b = math.sqrt(y)
    print(b, b * b)

  结果还是一样的：

29.79932885150268 888.0
29.79932885150268 888.0

  为什么呢？因为$\frac{1}{2}(x+\frac{y}{x})$，如果$x$初始化大了，$1/2$会缩小它，如果$x$小了，$y/x$会扩大它，最终接近于$\sqrt{x}$。但是初始值还是重要的，毕竟可以减少一点迭代次数，选择$x=y$这个初始值就非常不错。