【编程孔乙己】之开方倒数的几种写法孔乙己用指甲蘸了酒，问道："我便考你一考，回字有四样写法，你知道么？"不久前，美团的面

孔乙己用指甲蘸了酒，问道："我便考你一考，回字有四样写法，你知道么？ "

不久前，美团的面试官问我，"我便考你一考，你知道开平方么？便写一个罢 "
我心里只隐约记得大学课程的《数值分析》里介绍了一种叫做牛顿迭代的算法进行计算，算法细节却早已忘光。于是我便写了一个愧对职业素养的实现：暴力破解

1.暴力破解
思路很简单，我随便找一个初始值，求出平方与目标比对一下，大了就减小，小了就增大。反复尝试总会有对的时候

    /**
     * 暴力破解，傻瓜方案
     */
    @Test
    public void foolSqrt() {
        double target = 9;
        double point = target / 2;
        int times = 0;
        double precision = 1e-2;
        double gap = 1;
        while (Math.abs(gap) > precision) {
            if ((target - point * point) > 0) {
                point += precision;
            } else {
                point -= precision;
            }
            gap = target - point * point;
            times++;
        }
    }

为了求出9的平方根，尝试 : 151次 , point : 3.000000000000032 , point的平方 : 9.000000000000192 , point平方与target的差值 : -1.9184653865522705E-13

2.牛顿迭代

牛顿迭代其实不仅仅是用来做开平方的，而是快速逼近求得方程的逼近解。以下内容来自果壳文章《求牛顿开方法的算法及其原理，此算法能开任意次方吗?》

假设方程在附近有一个根，那么根据f(x)在x0附近的值和斜率，就能估计f(x)和x轴的交点用。迭代式子：依次计算、、、……，那么序列将无限逼近方程的根。

用牛顿迭代法开平方】令：

所以f(x)的一次导是：

牛顿迭代式：
随便一个迭代的初始值，例如，代入上面的式子迭代。例如计算，即a=2。

……
计算器上可给出

根据果壳文章以上阐述，由此我们可以写出如下代码：

    /**
     * 牛顿迭代开平方
     */
    @Test
    public void newtonSqrt() {
        int target = 9;
        int point = 1;
        double precision = 1e-2;
        int times = 0;
        while (Math.abs(target - point * point) > precision) {
            point = (point + target / point) / 2;
            times ++;
        }
    }

求9的平方根，尝试 : 3次 , point : 3 , point的平方 : 9 , point平方与target的差值 : 0

到此为此，似乎问题已经得到了圆满的答案。但只有一个正确答案如何对得起编程孔乙己的title。紧接着我就发现了更有趣的解法

3.魔法数开方

暴雪是一家神奇的公司，它的代码总是能给人惊喜，例如 One-way Hash 被业内赞为最优秀的hash表改良算法（timer33算法应该才是最快的，One-way Hash 在Timer33的基础上改用两次Hash计算来比对Key是否一致，在数据结构没有根本变化的情况下，理论上来说应该更慢）。

下面是暴雪的开平方算法源码，它对上述的牛顿迭代做了一些改进：

#include <math.h>  
float InvSqrt(float x)  
{  
 float xhalf = 0.5f*x;  
 int i = *(int*)&x; // get bits for floating VALUE  
 i = 0x5f375a86- (i>>1); // gives initial guess y0  
 x = *(float*)&i; // convert bits BACK to float  
 x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy  
 return x;  
}  
  
int main()  
{  
  printf("%lf ",1/InvSqrt(3));  
  
   return 0;  
}

这个算法的牛逼之处在于，经典的牛顿迭代里猜测值是随便取的一个值，而这段代码里却用了一个魔法数0x5f3759df来作为猜测值。这个算法只用了一次位运算，根本没有多次迭代就得到了结果，比直接用sqrt(n)快了大约四倍。

但是，我要说但是了，这个算法并不是暴雪原创的。它更早的时候出现在雷神之锤3 这个游戏的3D引擎代码里，作者是约翰-卡马克（John Carmack），江湖人称卡神。

float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;

x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}

注意看代码 i = 0x5f3759df - ( i >> 1 ); // what the fuck?很显然，卡神也不知道这个魔法数是哪来的。根据gamedev一群大神的挖掘，发现早在70年代NASA的代码里就有了这个魔法数。

这个魔法数 0x5f3759df 。它是从哪来的，又起到了什么作用？
普渡大学的数学家Chris Lomont在gamedev上看到这个讨论以后觉得有趣，决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。为此写下一篇论文， Fast Inverse Square Root 并通过暴力穷举的方式，得出一个更好的开方猜测值：0x5f375a86

为什么这个值是0x5f3759df，Chris Lomont的论文Fast Inverse Square Root给出了答案（点我看WIKI）知乎上有同学给出了汉化版答案（点我看知乎原文）：

单精度的浮点数结构是这样：

那么现在，每个正浮点数 y 可以用尾数和指数的形式写成 2^e(1+m)

，其中 m 是尾数部分，取值范围是 [0, 1)

；e 是指数部分，一个整数。每个浮点数所对应的「整数形式」则可以用 EL+M

表示，其中 L 是指数部分需要的位移次数（用 2 的幂表示），E 和 M 是指数部分和小数部分的整数版本。两者之间的关系是
$\left\{ \begin{array}{ll} E=e + B \\ M=m L \end{array} \right.$
对单精度浮点而言， $L=2^{23},B=127$ 。

考虑对数 $\log_2 y=e+\log_2(1+m)$ ，由于 $m\in\left[0, 1\right)$ ， $\log_2(1+m)\in\left[0,1\right)$ ，取近似 $\log_2(1+m)\approx m+\delta$ ，可以算出整体「偏差」最小的 $\delta=0.0430357$ ，此时两者基本相当。因此我们可以说

$\log_2 y\approx e+m+\delta$ ……………… (1)

那么，对于 y 的整数形式 Y 而言，展开并带代入 (1) 有：
$Y &=& amp; EL+M\\ &=& amp;L(e+B+m) \\ & \approx & L(\log_2 y + B - \delta)$
即

$\log_2 y\approx {Y \over L}-(B-\delta)$ ………………(2)

那么对于 $y '={1 \over \sqrt{y}}$ 来说， $\log_2{y'}\approx -{1\over 2}\log_2 y$ ，代入 (2) 得：
${Y ' \over L}-(B-\delta)={1 \over 2}\left[(B-\delta)-{Y\over L}\right]$
解得
$Y'={3 \over 2}L(B-\delta)-{1 \over 2}Y$ 。
这个就是代码



 i = 0x5f3759df - ( i >> 1 );

的秘密所在。

程序员多学点数学，总是会有用的，与君共勉附录：java版算法实现

public static double fastInvertSqrt(double number) {
    long i;
    double x2, y;
    final double threehalfs = 1.5F;
    x2 = number * 0.5F;
    y = number;
    // evil floating point bit level hacking
    i = Double.doubleToLongBits(y);
    // gives initial guess y0
    i = 0x5fe6eb50c7b537a9L - (i >> 1);
    y = Double.longBitsToDouble(i);
    // 1st iteration
    y = y * (threehalfs - (x2 * y * y));
    // 2nd iteration
    y = y * (threehalfs - (x2 * y * y));

    return 1 / y;
}

引用文献：

[1].果壳《求牛顿开方法的算法及其原理，此算法能开任意次方吗? 》www.guokr.com/question/46…

[2].维基百科《Fast inverse square root》en.wikipedia.org/wiki/Fast_i…

[3].知乎《0x5f3759df这个快速开方中的常数的数学依据是什么？》www.zhihu.com/question/20…