【转】开平方取倒数

Quake III 公开源码后,有人在game/code/q_math.c里发现了这样一段代码。它的作用是将一个数开平方并取倒,经测试这段代码比(float)(1.0/sqrt(x))快4倍。具体代码为:

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

函数返回1/sqrt(x),这个函数在图像处理中比sqrt(x)更有用。这个简洁的函数,最核心,也是最让人费解的,就是标注了“what the fuck?”的一句

i = 0x5f3759df( i >> 1 );

再加上

y  = y * ( threehalfs – ( x2 * y * y ) );

两句话就完成了开方运算!
算法的原理其实不复杂,就是牛顿迭代法,用x-f(x)/f’(x)来不断的逼近f(x)=a的根。一般的求平方根都是这么循环迭代算的,但是整个程序中选择了一个神秘的常数0x5f3759df 来计算那个猜测值,算出的值非常接近1/sqrt(n),这样我们只需要2次牛顿迭代就可以达到我们所需要的精度。
数学家Chris Lomont看了以后觉得有趣,决定要研究一下弄出来的这个猜测值有什么奥秘。在精心研究之后从理论上也推导出一个最佳猜测值, 0x5f37642f。Lomont计算出结果以后非常满意,于是拿自己计算出的起始值和神秘数字做比赛,看看谁的数字能够更快更精确的求得平方根。结果是Lomont输了…
via http://www.biaodianfu.com/0x5f3759df.html