平方根算法之一(翻译)

LiHengFeng posted @ 2008年3月19日 04:48 in 算法 , 4121 阅读

本blog所有文章，除声明了作者外，均为原创。欢迎转载，转载请注明作者。

平方根算法

原作者： Paul Hsieh

原文出处：http://www.azillionmonkeys.com/qed/sqroot.html

翻译：lhf

这里有关于平方根的一切，也许是你曾经想知道但却不好意思问的。

一、计算方法

大部分编程语言，如C和X86汇编都带有平方根函数。但是有许多原因让我们想要一个软件的模拟算法来计算平方根。其中一个原因是既不想在计算精度上妥协，又要更好的执行性能，另外的原因还有在某些系统上平方根运算没有得到很好的支持。

当遇到平方根问题时，我们可能首先会采用这样的估算方法：用y平方根的一个估计值x_est用作为起始值，然后通过下面的公式反复迭代：

x_est ← (x_est + y/x_est) / 2

公式背后的道理很简单。如果估计值x_est很小，那么y/x_est果将会很大，反之亦然。因此，我们可以期待： x_est < y < y/x_est（或者反过来x_est > y > y/x_est），并希望这两个数的平均值是一个更好的平方根估计。当0 < x_est < y成立时，确实是这样的（另一方面，如果x_{est，容易知道在有限次迭代后估计值会小于y+2, 因此在n次进一步迭代后会小于y/2 + 1/2ⁿ）。你可以在函数f(y, x) = y – x²上使用牛顿方法得到这个结论。} > y

现在，让我们取：x_est = (√y) * (1 + e)。在一次迭代后得到：x_est = (√y) * (1 + e²/(2*(1+e)))。对于一个比较小的e值，每次迭代误差由e变为e²/2，因此在每次迭代后，估计值的有效数字或者位数都会加倍。这就是说，一旦估计值的最高位是正确的，对53个二进制位的精度来说，仅需要不超过六次迭代。对任何e > -1的值来说，迭代值总是更加精确的，这意味着对任何非零的初始值，估计值总是收敛的。

这里要说的第二个方法是一种使用增量来估计平方根的方法。当需要手算平方根，或者只需要较低精度时，我们通常会采用这种方法。

假定我们有一个对y平方根的估计值：x_est，可以用下面的误差函数来测量它的精度：

err(x,y) = y - x_est²

现在考虑估计值x+p的误差：

err(x_est+p,y) = y – (x_est+p)²
= y – x² – 2*x_est*p – p²
= err(x_est,y) – p*(2* x_est + p)

我们希望err(x+p,y)能尽可能的接近0。如果假定x比p要大的多，则我们可以从下面的等式中求出p：

0 = err(x_est,y) – p*(2*x_est)

求出p = err(x_est,y)/(2*x_est)。因此，如果我们用一个合理的初始估计值x_est开始，就可以通过下面的步骤迭代产生一个越来越好的估计值。

1、x_est 是一个y平方根的估计值

2、e = y – x_est²

3、p = e/(2*x_est)

4、e = e – p*(2*x_est + p)

5、x_est = x_est + p

6、如果 x_est 不够好，回到第3步。

7、output x_est

取：x_est = (√y) * (1 + e). 在一次迭代后同样得到：x_est = (√y) * (1 + e²/(2*(1+e)))。因此，实际上这个方法的效果和上一个方法是完全相同的。这个方法的优点是将第一个方法重新表达成一种新的形式，在这种形式下，我们可以在第3步有些控制权。例如：如果第3步的除法只精确到一个数字或一位，算法的每一步将精确收敛到一个数字或一位的精度。因此，从除法运算花费的代价或者用指定的低精度除法得到近似值的精度来看，这个方法比第一个更好。

一个简单的问题是，能否在算法中完全不使用除法？答案很有意思。我们不能不用除法而直接求平方根，但却可以求平方根的倒数。如果求出了1/√y，那么 y * (1/√y) 就是我们想要的结果。我们从平方根倒数的估计值x_est开始，用下面的公式迭代：

x_est ← x_est*(3 - y*x_est²) / 2

这个公式可以在函数f(x) = y – 1/x²上应用牛顿方法得到。从公式中可以看到，如果估计值x_est太小，那么(3 - y*x_est²)将大于2。如果估计值太大，则(3 - y*x_est²)将小于2。但不幸的是，现在算法要成功的收敛，对初始估计值的依赖稍微有的紧密了些。具体来说，要求估计值满足：0 < y*x_est² < 3。

令 x_est = (√y) * (1 + e)。一次迭代后得到：x_est =(1/√y) * (1 – (3/2)*e²- e³/2)。对于比较小的e，每次迭代误差由e 变为-(3/2)*e²，因此每一轮迭代，有效的数字或者位数几乎会加倍。尽管需要更多的迭代次数，但完全不使用除法会很容易弥补这一点。

现在的问题是，如何得到一个可靠的估计值满足：0 < x_est< √3/√y，当然这里不能使用任何平方根函数。

很幸运的是，我们可以利用IEEE-754浮点数格式。任何非0，inf，-inf，NaN的数都可以表示为如下的格式：

x = sign * (1 + M) * 2^E, 其中0 ≤ M < 1, E为整数，sign 为 -1或者 1.

假定sign为1，则平方根可以大致表示为：

x = √(1 + M)* 2^(E/2)

使用(1+x)^1/2的级数展开，可以得到：

x ≈ (1 + M/2 + O(M²)) * 2^{floor (E/2)} * (1 + (√2-1) * Ind(E is odd))

我们知道，对一个小的x，有1/(1+x) ≈ 1-x+O(x²)，因此可以取倒数如下：

x ≈ (1 - M/2 + O(M²)) * 2 ^{-floor (E/2)}* (1 – (1-1/√2) * Ind(E is odd))

对上式，如果不关心2的幂部分，可以注意到它基本上是独立于指数部分的。我们可以检查一下估计值的精度，比如对1到4的范围内进行估计。在1到2-epsilon范围内，最大误差在2-epsilon处，公式产生的估计值是0.5，真实值是0.707。在2到4-epsilon的范围内，最大误差在4-epsilon处，此时估计值是0.353，真实值是0.5。这样就很容易得到满足0 < x_est< √3 / √y的估计值。

现在，我们剩余的难题是，初始估计值可能会太好，我们也许想要简化它通过避免乘以(1 – (1-1/√2) * Ind(E is odd))的因子。我们注意到上面的公式可以改写为如下：

x ≈ (1 - M/2 + 0.5*Ind(E is even)) * 2 ^{-floor (E/2)}

我们估计的范围仍然有效。那么为什么会想到将公式改写为这样的形式呢？因为IEEE-754格式是这样的：

V = [sign (=0)] [biased-exponent (=Bias + E)] [mantissa (= M)]

考虑到作为对恰当的位做一个类似整数的操作，我们希望有一个常数K，使得K-((unsigned long)V)>>1产生一个IEEE-745格式的浮点数，恰好对应于我们要找的估计值。注意到尾数和指数都划分为两部分，如果基于偏离的指数部分是奇数，会附加一个0.5到尾数上。因子K仍然需要正确的指数基准以便产生正确的表示格式。实际上，对IEEE-754的64位格式，我们发现这个常数是0xBFCD460000000000，对32的是0xBE6F0000。可以看到IEEE-754的设计规格以前事先考虑了可以很好的应用这些技巧。