求平方根算法 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\)。算法描述如下：
- 初始化 \(x = 1\);
- 计算 \(\Delta_k = |x_k – x_{k-1}|\);
- 若\(\Delta_k < \epsilon\)，返回\(x\)，算法结束;
- \(x \leftarrow \frac{x}{2} + \frac{m}{2 x}\);
- 返回第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]\)区间：
- 初始化 \(u=0\);
- 若 \(0.5 \leq m \leq 2\),返回 \(u\)，算法结束；
- 若 \(m > 2\): \(m \leftarrow m / 4\), \(u \leftarrow u+1\);
- 若 \(m < 1\): \(m \leftarrow m * 4\), \(u \leftarrow u-1\);
- 返回第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的迭代公式了。