开方的几种写法

1.牛顿迭代

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

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

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

    /**
     * 牛顿迭代开平方
     */

public static double newtonInvertSqrt(double target) {
    double point = 1;
    double precision = 1e-2;
    while (Math.abs(target - point * point) > precision) {
        point = (point + target / point) / 2;    

    }

    return point;

}

到此为此，似乎问题已经得到了圆满的答案，但紧接着我就发现了更有趣的解法

2.魔法数开方倒数

暴雪是一家神奇的公司，它的代码总是能给人惊喜，例如 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）知乎上有同学给出了汉化版答案（点我看知乎原文）：

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

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;
}