【转】代码优化之-优化开根（卡马克算法）

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

          

        由这个公式我们就很清楚地明白代码y=y*(threehalfs-(x2*y*y))的含义，这其实就是执行了单次牛顿迭代。为啥只执行了单次迭代就完事了呢？因为单次迭代的精度已经达到相当高的程度。 

         为什么单次迭代就可以达到精度要求呢？根据之前的分析我们可以知道，最根本的原因就是选择的初值非常接近精确解。而估计初始解的关键就是下面这句代码：

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

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

          float型浮点数和对应的int型整数之间的关系给出一个公式：

         

     其中，Ix表示float型浮点数地址强转后的int型整数，L=2^23，x是原始的浮点数（尚未表示成float类型），B=127，是一个无穷小量。化简一下上述公式我们得到：

 

            

        有了这个公式我们就可以推导初始解的由来了。要求，我们可以将其等价转化成，然后代入上面的公式我们就得到：

    这个公式就是神奇操作的数学表示，公式中只有是未知量，其它都已知。的值没有好的求解方法，数学家通过暴力搜索加实验的方法求得最优值为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

y₁ =y(1.5-2.5625y^2)≈ 0.441593890

 

 

参考来源：夏风习习的博客->卡马克算法，http://blog.163.com/lxd007_2005/blog/static/405618252015112410210140/

 

 

 

 

附：整数开根处理函数

 

/**************************************************************
 * 函数名：Q_rsqrt
 * 描  述：开根号函数（整型）
 * 输  入：uint32_t radicand，被开方数
 * 输  出：uint32_t result_sqrt，开方结果
 * 说  明：返回开方后的整数结果，运行16次。
**************************************************************/
uint32_t Q_rsqrt(uint32_t radicand)
{
    uint32_t i;    
    uint32_t result_sqrt;            // 开方结果
    uint32_t radicand_top,radicand_rmng;    
    if(radicand == 0)
    {
        return 0;
    }
    result_sqrt = 0;
    radicand_top = (radicand >> 30);
    radicand <<= 2;
    if(radicand_top > 1)
    {
        result_sqrt++;
        radicand_top -= result_sqrt;
    }
    for(i=15; i>0; i--)
    {
        result_sqrt <<= 1;
        radicand_top <<= 2;
        radicand_top += (radicand >> 30);
        radicand_rmng = result_sqrt;
        radicand_rmng = (radicand_rmng << 1) + 1;
        radicand <<= 2;
        if(radicand_top >= radicand_rmng)
        {
            radicand_top -= radicand_rmng;
            result_sqrt++;
        }
    }    
    return result_sqrt;
}