【概述】
越底层的函数,调用越频繁,那么最底层的数学运算函数的优化至关重要。
当求逆平方根时,一般做法都是用函数返回 1/sqrt(x),但在 雷神之锤3 中,有一快速求逆平方根的算法。
【知识储备】
1.单精度浮点数的存储
在计算机中,单精度浮点数使用 32 位来存储,其中,最高位为符号位,后面 8 位为整数位 ,代表浮点数的指数,再后面 23 位代表小数部分
,依次表示
、
、...
因此,如果 x 是一个正浮点数,则有:
如果想将一浮点数转为整数形式,则需要做如下运算:
其中,L 是指数部分唯一需要的次数,M 是小数部分对应的整数版本,B 为 127
2.牛顿迭代法
牛顿迭代法,是求解任意连续函数的根值的一种方法。
对于如上函数,现要求这个函数的根:首先猜一个 ,假设它是函数的解,但由于其实际不是,因此需要将这个解迭代,使其逼近真正的解。
在 处作其切线,求得切线方程:
,并求切线的根,可以发现对于真正的解已经逼近了一步。
推广到 n,继续迭代,就足够逼近真正的解:
此时, 可被一个统一的函数来表示:
令 ε 为当前解与真正解 r 的距离,则:
综合方程,可得:
因此,只要 ε 小于某个特定的值,则可认为此时的 与方程的解十分接近
【源码】
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;
}
【分析】
如果要求一个浮点数的平方根倒数,一般求法为:
将其转化为关于 y 的方程,有:
转换为牛顿迭代法的方程,有:
此时,对原方程两边同取 2 的对数,有:
由于 ,则在这个区间内,可以近似取:
根据方差的计算,当 σ=0.0430357 时,整体的偏差是最小的,此时上面的等号两边应该相当。
将上述公式整合,最终的 可以写成:
则:
其中:
最终,写成代码就是:
i = 0x5f3759df - ( i >> 1 );
【关于源码】
求逆平方根的算法来自 雷神之锤3(quake3) 的作者卡马克,该算法并不复杂,其核心就是用牛顿迭代法来不断逼近,但卡马克真正厉害的地方,在于他选择了一个十分神秘的常数:0x5f3759df 来计算猜测值,于是第一次牛顿迭代算出的值已经非常接近 ,这样仅需两次牛顿迭代就可达到所需精度。
普渡大学的数学家 Chris Lomont 看了以后觉得有趣,决定要研究卡马克弄出来的这个猜测值有什么奥秘。
Lomont 在精心研究之后从理论上也推导出一个最佳猜测值:0x5f37642f,和卡马克的数字非常接近。
传奇并没有在这里结束。
Lomont 计算出结果以后非常满意,于是拿自己计算出的起始值和卡马克的神秘数字做比较,看哪个数字能够更快更精确的求得逆平方根,结果是卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。
最后 Lomont 采用暴力方法一个数字一个数字试过来,终于找到一个比卡马克数字要好上那么一丁点的数字,虽然实际上这两个数字所产生的结果非常近似,这个暴力得出的数字是:0x5f375a86。
Lomont 为此写下一篇论文,"Fast Inverse Square Root"。
论文下载地址:
http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
http://www.matrix67.com/data/InvSqrt.pdf