神奇的0x5f3759df

会数学的程序员都是艺术家，第一篇文章致敬John Carmack的wtf magic number 0x5f3759df  —— 题记

　　以前上网偶遇QuakeIII的完全平方根倒数的代码，就被这magic number后的注释吸引了，以致于很长一段时间都喜欢在code review时用0x5f3759df作为我wtf的文明嘴替来comment。直到有一天突然好奇这玩意到底是啥，就翻了一下源码去分析，折腾了1个多小时才搞明白，此后除了对作者John Carmack的敬佩，就只剩想被0x5f3759df comment？你不配

 

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

 

　　代码见上方block，分享下到时思路。正着看下来估计谁看谁懵逼，那就反着吧，13行可以确认y就是计算结果res，再看11、12行一样的执行，等号两边都有res，再加上注释里的iteration，牛顿迭代法没跑了，手推一波

 

　　完美，对上了。那就这堆莫名其妙的7-10行的意图清晰了，就是搞个初值用来估算。下面思路有了，查下float的存储方式，然后按着0x5f3759df前后代码手推一波，应该就能知道wtf是wtf了

 

　　bingo，其实就是将log(1+x)在x∈(0,1)时近似成了x+0.0450465679169。当然，这有俩限制一是不能搅屎棍输入负数，二是输入的数E不能为0即不能小于2^-126。

　　第二个限制是不是真实不行呢，懒得推了，直接编程随机测试下

#include <stdio.h>
#include <stdlib.h>

#include <math.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?
    y = * ( float * ) &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()
{
    long tmpFloatMaker = (0U << 31) | (0U << 23) | (0x1 << 0);
    float num = *(float*)&tmpFloatMaker;

    printf("num = %.60f\n",num);
    printf("Q_rsqrt(num) = %.20f\n",Q_rsqrt(num));
    printf("1/sqrt(num) = %.20f\n",1/sqrt(num));

    return 0;
}

 

　　显然，输入的数不能小于2^-126，因为E全0后float对M的解析不再是+1.0而是+0.0，2的幂也变为-126而不是预期的-127。

　　好了，现在正着看回来，这串代码的思路就是先用一个很神奇近似得到一个近似解，然后采用(1/*+1*/)次牛顿迭代，得到更精确点的近似解。至于那个神奇近似方法，正着真不知道怎么想的，颇有拉马努金的感觉，0x5f3759df