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;

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

这是在《雷神之锤 III 竞技场》源代码中的一个函数，功能是计算 平方根倒数1x\dfrac{1}{\sqrt{x}}\\.

在 3D 游戏场景的计算机图形学中，要求取照明和投影的光照与反射效果，就需要对一个曲面上的多个点坐标计算平方根倒数，即

1x2+y2+z2\dfrac{1}{\sqrt{x^2 + y^2 + z^2}}

浮点数的表示·

浮点数在线转换工具：IEEE-754 Floating Point Converter

在 C 语言中，浮点数表示采用 IEEE 754 标准。对一个 32 位二进制浮点数，从高到低分成三部分

  • [31:31]：表示符号位(Sign, S)；
  • [30:23]：表示指数(Exponent, E)；
  • [22:0]：表示尾数(Mantissa, M)；

以小数 3.143.14 为例，

  • 非负数，符号位为 S=0S=0；
  • 8 比特的指数位表示 [0,255][0, 255] 区间内的整数，为了表示小于 11 的浮点数，指数需要取负数，所以将取值修改在 [−127,128][-127, 128] 这个区间；因为 21<3.14<222^1 < 3.14 < 2^2，所以 E=128E=128；
  • 23 比特偏移量，取值范围在区间 [0,223][0, 2^{23}]，将 [2n,2n+1][2^n, 2^{n+1}] 这个区间分成了 2232^{23} 段，MM 表示从 2n2^n 开始的第 MM 段，M=3.14−24−2=0.57,0.57×223=4781507M = \dfrac{3.14-2}{4-2}=0.57, 0.57\times 2^{23} = 4781507，

最终得到 32 位二进制转换为浮点数的计算公式为

(−1)S×(1+M223)×2(E−127)(-1)^S \times (1 + \dfrac{M}{2^{23}}) \times 2^{(E - 127)}

其中 S=0,M=4781507,E=128S=0, M=4781507, E=128，得到 3.143.14 的二进制表示为

0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 1 1 0 1 0 1 1 1 0 0 0 0 1 1

带入公式得到近似值为 3.14000010490417480468753.1400001049041748046875

平方根倒数计算·

1x\dfrac{1}{\sqrt{x}}\\的定义域与值域均为[0,+∞)[0, +\infty)，S=0S=0，于是浮点数的计算公式可以简化为

(1+M223)×2(E−127)(1 + \dfrac{M}{2^{23}}) \times 2^{(E - 127)}

该 32 比特二进制数表示的整型为 M+223EM + 2^{23}E。

y=1xy = \dfrac{1}{\sqrt{x}}

两边同时取对数，得到

log⁡2y=−12log⁡2x\log_2{y} = -\dfrac{1}{2}\log_2{x}

将 x,yx, y 替换为 (1+M223)×2(E−127)(1 + \dfrac{M}{2^{23}}) \times 2^{(E - 127)}，

log⁡2(1+My223)+Ey−127=−12(log⁡2(1+Mx223)+Ex−127)\log_2{(1+\dfrac{M_y}{2^{23}})} + E_y-127 = -\dfrac{1}{2}(\log_2{(1+\dfrac{M_x}{2^{23}})} + E_x-127)

令 log⁡2(1+M223)≈M223+σ\log_2{(1 + \dfrac{M}{2^{23}})} \approx \dfrac{M}{2^{23}} + \sigma，则

My223+σ+Ey−127≈−12(Mx223+σ+Ex−127)\dfrac{M_y}{2^{23}} + \sigma + E_y-127 \approx -\dfrac{1}{2}(\dfrac{M_x}{2^{23}} + \sigma + E_x-127)

整理得

My223+Ey≈−12(Mx223+Ex)−32(σ−127)My+223Ey≈3×2232(127−σ)−12(Mx+223Ex)\begin{aligned} \dfrac{M_y}{2^{23}} + E_y &\approx -\dfrac{1}{2}(\dfrac{M_x}{2^{23}} + E_x) -\dfrac{3}{2}(\sigma-127)\\ M_y + 2^{23}E_y &\approx \dfrac{3\times 2^{23}}{2}(127 - \sigma) - \dfrac{1}{2}(M_x + 2^{23}E_x)\\ \end{aligned}

其中 My+223EyM_y + 2^{23}E_y 为表示浮点数 yy 的 32 比特二进制的整型，Mx+223ExM_x + 2^{23}E_x 为表示浮点数 xx 的 32 比特二进制的整型，3×2232(127−σ)\dfrac{3\times 2^{23}}{2}(127 - \sigma) 是一个常数。

文首程序中 i = * ( long * ) &y; 首先将浮点数的 32 比特二进制转换为整型，然后 i = 0x5f3759df - ( i >> 1 ) 与刚刚推导出的公式对应，3×2232(127−σ)=0x5f3759df\dfrac{3\times 2^{23}}{2}(127 - \sigma)=\text{0x5f3759df}。

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;

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

对于一个输入的浮点数 xx，首先得到首次近似值 0x5f3759df−12x\text{0x5f3759df} - \dfrac{1}{2}x，然后代入牛顿迭代公式 xn+1=xn−f(xn)f′(xn)x_{n+1}=x_n - \dfrac{f(x_n)}{f'(x_n)}。

对于 y=1xy=\dfrac{1}{\sqrt{x}}，f(y)=1y2−xf(y)=\dfrac{1}{y^2}-x，则

yn+1=yn−1yn2−x−2yn3=yn(1.5−0.5xyn2)y_{n+1} = y_n - \dfrac{\frac{1}{y_n^2}-x}{-\frac{2}{y_n^3}} = y_n(1.5-0.5xy_n^2)

对应代码 y = y * ( threehalfs - ( x2 * y * y ) );.

只搞懂了这段代码的数学原理，还是不知道 0x5f3759df\text{0x5f3759df} 这个数从何而来……

参考资料·