开方的几种写法

1.牛顿迭代

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

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

用牛顿迭代法开平方】令：
 
所以f(x)的一次导是：

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



……
计算器上可给出

 

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

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

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;

}

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

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

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）知乎上有同学给出了汉化版答案（点我看知乎原文）：

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

那么现在，每个正浮点数 y 可以用尾数和指数的形式写成 ，其中 m 是尾数部分，取值范围是 ；e 是指数部分，一个整数。每个浮点数所对应的「整数形式」则可以用表示，其中 L 是指数部分需要的位移次数（用 2 的幂表示），E 和 M 是指数部分和小数部分的整数版本。两者之间的关系是

对单精度浮点而言，。

考虑对数 ，由于，，取近似，可以算出整体「偏差」最小的，此时两者基本相当。因此我们可以说 

……………… (1)

那么，对于 y 的整数形式 Y 而言，展开并带代入 (1) 有：

即

………………(2)

那么对于 来说，，代入 (2) 得：

解得
。

这个就是代码

 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].果壳《求牛顿开方法的算法及其原理，此算法能开任意次方吗? 》http://www.guokr.com/question/461510/

[2].维基百科《Fast inverse square root》https://en.wikipedia.org/wiki/Fast_inverse_square_root

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