在电影《微微一笑很倾城》中,肖奈大神在玻璃上写了一堆公式,提到平方根倒数速算算法,这个到底是一个什么算法?笔者看电影的时候打开手机学了一下,发现该算法的作者真乃神人!今天有空,就把该算法写一写。
在3D图形编程中,经常要求平方根的倒数,即1/Sqrt(x),如果用一般的代码(float)(1.0/sqrt(x)),,精度高,但是非常慢;我们需要一个快速,而又足够高精度的算法;著名游戏《雷神之锤III》的代码在2002年左右被披露,人们发现了一段用于快速计算平方根倒数的代码,下面是整理后的代码(去掉了一些宏定义)。
float InvSqrt (float x){float xhalf = 0.5f*x;int i = *(int*)&x; // get bits for floating valuei = 0x5f3759df - (i >> 1);// L1:gives initial guess y0x = *(float*)&i; //l2:convert bits back to floatx = x*(1.5f - xhalf*x*x); //l3:Newton step, repeating increases accuracyreturn x;}
该程序运行效率极高,经测试,基本是使用直接开根号求倒数程序的4倍速度!一时间惊为天人。那么这段代码到底怎么理解?为什么中间出现了0x5f3759df这样一个完全无厘头的magic number?
--------------------------------------------------------------------------------------------------------------------------------------------------
该算法的本质其实就是牛顿迭代法(Newton-Raphson Method,简称NR)。NR是一种求方程的近似根的方法。首先要估计一个与方程的根比较靠近的数值,然后根据公式推算下一个更加近似的数值,不断重复直到可以获得满意的精度。其公式如下:
函数:y=f(x)
其一阶导数为:y'=f'(x)
则方程:f(x)=0 的第n+1个近似根为
x[n+1] = x[n] - f(x[n]) / f'(x[n])
--------------------------------------------------------------------------------------------------------------------------------------------------
求一个数a的平方根的倒数,实际就是求方程f(x)=1/(x^2)-a=0的解;将该方程按牛顿迭代法的公式展开为:
x[n+1]=x[n]*(3/2-a/2*x[n]*x[n])
也就是上面代码蓝色部分的第L3行;所以很明显,这段代码就是进行一次迭代的牛顿迭代法。只进行一次迭代就结束程序,又要保证精度,也就是说,牛顿迭代法的初始值要非常精准,才能在一次迭代后完成计算。整个代码最精彩的就是L1行了,i = 0x5f3759df - (i >> 1);到底是什么意思?为什么初始值可以这样算呢?了解该代码,需要先了解一些float point和fix point的表示方法[2]:
一个浮点数float是由32位二进制位表示的有理数,分为三部分。其中符号占1位,表示正负,记为Si;指数占接下来的8位,表示经过偏移处理后的指数,即实际表示E(如图中为124),需要偏移B(图中为2的8次方减1,127。B为一个固定值),最后得指数值为E-B;有效数字(除最高位1以外,因为前面有指数,所以一个数肯定可以表示成1.xxxxx * 2^(kkk),只保留xxxxx)占剩下的23位,记为m(0<m<1),图中的。
所以浮点数的结构公式为:, 图中
整数fix point的表示相对简单,符号占1位,数值占剩下的31位。如果用上图的浮点数字节序列来表示整数,那么,即;平方根倒数函数仅能处理正数,所以符号位均为0。
对于同样的32位二进制数码,若为浮点数表示时实际数值为,而若为整数表示时实际数值则为,其中,这里n=32,b=8。式子中引入的新变量为:
------------------------------------等式1
,其中------------等式2
理解浮点数和整数的表示后,下面开始推导。
x的平方根倒数方程为:
两边取对数有:
因为浮点数可表示为:,所以也有,代入上式有:
由于,有,则在此可定义与x的关系为,表示误差,所以将代入上式得:
将等式1,等式2代入到上述方程中,有:
移项整理得:
又因为浮点规格存储的正浮点数x,若将其作为整数表示,则示值为:,所以x的平方根倒数的首次近似值的整数表示值为:
,
其中,,,n=32,b=8。
这个式子对应着源代码中的这一行:i = 0x5f3759df - ( i >> 1 );,然后将整数表示值换回表示浮点数:x = * ( float * ) &i;。这样就得到了浮点数的平方根倒数的近似值。
0x5f3759df 对应着R,即3/2(B-)L.当R为0x5f3759df时,有.
"现在不仅该算法的原作者不明,人们也仍无法明确当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验:他编写了一个函数,以在一个范围内遍历选取R值的方式将逼近误差降到最小,以此方法他计算出了线性近似的最优R值0x5f37642f(与代码中使用的0x5f3759df相当接近),但以之代入算法计算并进行一次牛顿迭代后,所得近似值与代入0x5f3759df的结果相比精度却仍略微更低。……在Charles McEniry的论文中,他使用了一种类似Lomont但更复杂的方法来优化R值:他最开始使用穷举搜索,所得结果与Lomont相同;而后他尝试用带权二分法寻找最优值,所得结果恰是代码中所使用的魔术数字0x5f3759df"---维基百科
参考资料:
[1]CHRIS LOMONT,FAST INVERSE SQUARE ROOT
[2]/heloowird/article/details/21862251
[3]/uid-9255716-id-107951.html
[4]/xiaoguohaha/article/details/21652643