回顾《雷神之锤III》平方根倒数算法

发表于 更新于 分类于

关于平方根倒数快速算法

《雷神之锤III》这款游戏的源码公布后，有人在里面发现了一个好玩的算法，是关于平方根倒数的算法。虽然该算法如今已经过时了，但是其中涉及的浮点数标准和牛顿迭代法还是挺值得一看的。

该算法一般用于求向量模长等地方，可以说在游戏建模中涉及物理还原这块用到这个平方根倒数还挺多的。

如果给我一个公式如下：

$$
f=\frac{1}{\sqrt{x}}
$$
我想我可能会写下

1
2
# include <math.h>
float y = sqrt(x);

但是在《雷神之锤III》里的算法并不是math.h里的算法。

源码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
float Q_rsqrt( float number )
{
	long i;
	float x2,y;
	const float threehalfs = 1.5F;
    
	×2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;				// evil floating point bit hack
	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, can be removed
	return y;
}

这个Q_rsqrt函数可以说是越看越怪，其中的0x5f3759df是哪里来的简直一头雾水。

IEEE754

首先需要了解浮点数的存储，IEEE二进制浮点数算术标准（IEEE 754），以下为例子，具体请自行了解。

其中第一位0正1负。

第一步：

1
i  = ( long )y;

将y强制转换为long格式是为了让y能够进行位操作，float并不能进行位操作。 （long 32bit）

但是转换会有误差产生，比如y=3.33在经过强制转换为long格式后就会变为3。同时内存小数部分也被打乱。

1
2
3.33	0	10000000	10101010001111010111000			//IEEE 754 
3		00000000	00000000	0000000	00000011

此时并不能达到要求，即对float格式的y进行位操作，毕竟内存中的位都改变了。

于是就有了

1
i  = * ( long * ) &y;				// evil floating point bit hack

这里首先通过&y获得了y的地址，然后通过(long *) 将改地址转换为long格式的地址，最后通过*获取地址内的值返回给i。

此时

上面这句话就是对第一步的理解，通俗来说就是原来IEEE 754的y的地址值是float格式，但是要保留地址内值不变的情况下告诉CPU这是long格式。

第二步：

1
i  = 0x5f3759df - ( i >>1 );		// what the fuck?

第三步：

牛顿迭代法
$$
f(x)=\frac{1}{x^{2}}-x
$$

$$
f^{\prime}(x)=\left.f^{\prime}(y)\right|_{y=x}=-\frac{2}{y^{3}}=-\frac{2}{x^{3}}
$$

$$
x_{\text {new }}=x-\frac{f(x)}{f^{\prime}(x)}=x-\frac{\frac{1}{x^{2}}-x}{-\frac{2}{x^{3}}}=\frac{3}{2} x-\frac{1}{2} x^{4}
$$