卡马克算法

原文地址：http://blog.163.com/lxd007_2005/blog/static/405618252015112410210140/

最原始的版本不是求开方，而是求开方倒数，也即。为啥这样，原因有二。首先，开方倒数在实际应用中比开方更常见，例如在游戏中经常会执行向量的归一化操作，而该操作就需要用到开方倒数。另一个原因就是开方倒数的牛顿迭代没有除法操作，因而会比先前的牛顿迭代（从Xi-1=1开始迭代）开方要快。

由这个公式我们就很清楚地明白代码y=y*(threehalfs-(x2*y*y))的含义，这其实就是执行了单次牛顿迭代。为啥只执行了单次迭代就完事了呢？因为单次迭代的精度已经达到相当高的程度。
为什么单次迭代就可以达到精度要求呢？根据之前的分析我们可以知道，最根本的原因就是选择的初值非常接近精确解。而估计初始解的关键就是下面这句代码：

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

正是由于这句代码，特别是其中的“magic number”使算法的初始解非常接近精确解。具体的原理是地址强转：首先将float类型的数直接进行地址转换转成int型（代码中long在32位机器上等价于int），然后对int型的值进行一个神奇的操作，最后再进行地址转换转成float类型就是很精确的初始解。
float型浮点数和对应的int型整数之间的关系给出一个公式

有了这个公式我们就可以推导初始解的由来了。要求，我们可以将其等价转化成，然后代入上面的公式我们就得到：
这个公式就是神奇操作的数学表示，公式中只有是未知量，其它都已知。的值没有好的求解方法，数学家通过暴力搜索加实验的方法求得最优值为0.0450466，此时第一项就对应0x5f3759df。但是后来经过更仔细的实验，大家发现用0x5f375a86可以获得更好的精度，所以后来就改用此数。
算法的最终目的是要对浮点数开平方，该算法性能非常高，而且精度也很高，三次迭代精度就和系统函数一样，但是速度只有系统函数sqrtf的十分之一不到，相当了得。

#include "stdio.h"

#include "conio.h"

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? 0x5f3759df or 0x5f375a86 什么该死的？ 卡马克算法 - dong - 北风寒*/

    y  = * ( float * ) &i;  /* 取长整型数i的地址，将其存储单元转换成浮点型，然后再把转换后的数取出来*/

    y  = y * ( threehalfs - ( x2 * y * y ) );   /* 1st iteration 第一次迭代*/   

    y  = y * ( threehalfs - ( x2 * y * y ) );   /*  2nd iteration, this can be removed 第二次迭代，能够移除*/ 

    return y;  

}  

int main()

{

    float n,z=1.0;

    printf("请输入一个需要求其平方根的数：");

    scanf("%f",&n);

    z=Q_rsqrt(n);

    printf("平方根为%f\n",1.0/z);

    getch(); 

    return 0;

}

举例：X=2^e(1+f)=5.125=2^2(1+0.28125)
Ix=EL+F=L(e+B+f)=2^23（2+127+0.28125）=2^23*10000001.01001=0(符号)10000001(阶码)
01001000000000000000000(尾数)（8388608*129.28125=1084489728）
Ix表示浮点数的整数表示，E=e+B表示IEEE阶码值，L=表示阶码的起始位置，F=Lf表示尾数的整数表示
=

12582912*(127-0.0450466)-1/2*1084489728=1597463007-542244864=1055218143=01111101 11001010101100111011111
y=*(float*)&i=2^(-2)*1.791805148124694824≈ 0.447951287
y1 =y(1.5-2.5625y^2)≈ 0.441593890