参考文章:一个Sqrt函数引发的血案
InvSqrt(value)函数相当于1.0/sqrt(value),只是一个简单地开平方取到数而已,但是下面这个代码竟然比1.0/sqrt(value)快了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;
}
简化后如下:
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
我们现在来一行一行地读一下代码。几年没碰过c了,赶紧复习了一下c。
首先第一行:
float xhalf = 0.5f*x;
这个就是把x乘0.5。
第二行来了:
int i = *(int*)&x; // get bits for floating VALUE
这一番操作差点把我难住了,毕竟我是搞python的,几乎不接触指针。
- &x
- 这是取x的地址
- (int*)&x
- 这是把x的地址转成int类型的指针。x本来是float类型的,在内存中存储自然是以float格式存储。但是经过这一操作后,就相当于告诉大家,这里存储的东西不归float管了,归我int管,以后你们要想从这读数据,必须用int格式读。
- 这好比你有一个png图片,你把后缀改成txt,然后双击打开,自然是用记事本打开,而且内容是乱码的。
- *(int*)&x
- 转换后在取出里面的内容,这时是以int类型取出。
我们一般不会把png改成txt,除非有某些东西不想让别人看到。
但是这里为什么有这个骚操作?(其实后面还有更骚的操作)
我不懂c,不知道这有什么意义,干脆直接手动给它算一遍。
令x=3,假设存储方式是32位的。
那么第一位是符号位,第2-9位是指数位,后面的是小数位。
3的二进制是11,写成指数是:
1.1
×
2
1
1.1\times2^1
1.1×21。
所以:
- 符号位是0。(0表示正,1表示负)
- 指数是1,这里因为指数有正有负,所以1+127=128,二进制为10000000。
- 小数位是
.
1
.1
.1(前面的1省略了,因为它总是1),即
2
−
1
2^{-1}
2−1,为10000000000000000000000(共23位)
- 顺便多解释一下,假如小数位是 . 1001101 .1001101 .1001101,即 2 − 1 + 2 − 4 + 2 − 5 + 2 − 7 + 2^{-1}+2^{-4}+2^{-5}+2^{-7}+ 2−1+2−4+2−5+2−7+,则为10011010000000000000000
所float类型的3以在内存中是:
01000000010000000000000000000000
01000000010000000000000000000000
01000000010000000000000000000000
现在,int霸占了它,于是我们用int格式来读它:
首位为符号位,0为正。
其他位为数字位,1000000010000000000000000000000的十进制是:
所以此时i=1077936128,我验证了一下,的确如此。
好了,我算出i了,但是不知道有什么意义,卡在这里没办法继续了。如果有朋友知道的话麻烦指点一下。
还是继续硬着头皮读下去吧:
更骚的第三行
i = 0x5f375a86- (i>>1); // gives initial guess y0
据说是用的牛顿法,这里是猜一个靠谱的y0
- i>>1
- i右移一位
- 即二进制为00100000001000000000000000000000
- 十进制为538968064(类似于除2再取整)
- 0x5f375a86
- 神秘的十六进制!
- 十进制:1597463174
- 二进制:1011111001101110101101010000110
。。。写不下去了