单位向量时需要用到平方根倒数,而计算单位向量在游戏引擎中会大量使用,属于底层代码,因此其效率将会直接影响游戏体验.
雷神之锤3中使用了以下代码
1 2 3 4 5 6 7 8 9 10 11 12 float Q_rsqrt (float number) long i; float x2, y; const float threehalfs = 1.5F ; x2 = number * 0.5F ; y = number; i = *(long *) &y; i = 0x5F3759DF - (i >> 1 ); y = *(float *) &i; y = y * (threehalfs - (x2 * y * y)); return y; }
本文将分析上述代码的含义
float类型总共占用32位,其内存结构采用科学计数法表示.第一位表示符号,接下来8位表示指数,后23位表示尾数
以4.25为例,其内存中的结构为0 10000001 00010000000000000000000
接下来的8位表示指数,其指在0-255之间,但是这样就无法表示负指数了,因此规定正指数第一位是1,负指数第一位是0,将这8位转换成10进制后减去127就是实际的指数.这里的10000001转换成10进制后是129,因此表示2 2 = 4 2^2 = 4 2 2 = 4
接下来的23位表示尾数.前面提到float在内存中是以科学计数法形式表示,在十进制中,科学计数法的个位数一定在1-9之间,因此在二进制中,科学计数法的个位数一定是1,既然一定是1,那就没有必要保存了,因此尾数0001实际上是1.0001,转换成十进制就是1.0625.
所以最终十进制数字就是1.0625 ⋅ 2 2 = 4.25 1.0625 \cdot 2^2 = 4.25 1 . 0 6 2 5 ⋅ 2 2 = 4 . 2 5
我们令8位指数转换成为十进制后为E,23位尾数转换成十进制后为M.
则原float数字可以表示为
( 1 + M 2 23 ) × 2 E − 127 \left(1+\frac{M}{2^{23}}\right) \times 2^{E-127} ( 1 + 2 2 3 M ) × 2 E − 1 2 7
取2为底的对数,得
log 2 ( 1 + M 2 23 ) + E − 127 \log_2 \left(1+\frac{M}{2^{23}}\right)+E-127 log 2 ( 1 + 2 2 3 M ) + E − 1 2 7
在上面的分析中,我们知道M 2 23 \LARGE \frac{M}{2^{23}} 2 2 3 M
M 2 23 + E + μ − 127 \frac{M}{2^{23}}+E+\mu-127 2 2 3 M + E + μ − 1 2 7
提出1 2 23 \frac{1}{2^{23}} 2 2 3 1
1 2 23 ( M + 2 23 ⋅ E ) + μ − 127 \frac{1}{2^{23}}\left(M+2^{23} \cdot E\right)+\mu-127 2 2 3 1 ( M + 2 2 3 ⋅ E ) + μ − 1 2 7
float无法进行位操作,而long可以,并且都是4字节,因此可以把float*转换成long*来进行位操作.
1 2 float y = number;long i = *(long *) &y;
现在开始正式计算y y y
则y y y
S = 1 y → log 2 S = − 1 2 log 2 y S=\frac{1}{\sqrt y}\rightarrow\log_2{S}=-\frac12\log_2y S = y 1 → log 2 S = − 2 1 log 2 y
因为上面的代码将y y y l o n g long l o n g M y M_y M y E y E_y E y M , E M,E M , E y y y y = 2 23 ⋅ E y + M y y=\textcolor{red}{2^{23}\cdot E_y+M_y} y = 2 2 3 ⋅ E y + M y S = 2 23 ⋅ E S + M S S=\textcolor{green}{2^{23}\cdot E_S+M_S} S = 2 2 3 ⋅ E S + M S
log 2 ( 2 23 ⋅ E S + M S ) = − 1 2 log 2 ( 2 23 ⋅ E y + M y ) \log_2{\left(\textcolor{green}{2^{23}\cdot E_S+M_S}\right)}=-\frac12\log_2{\left(\textcolor{red}{2^{23}\cdot E_y+M_y}\right)} log 2 ( 2 2 3 ⋅ E S + M S ) = − 2 1 log 2 ( 2 2 3 ⋅ E y + M y )
使用之前提到的近似替换方法,替换两个对数,得
1 2 23 ( 2 23 ⋅ E S + M S ) + μ − 127 = − 1 2 [ 1 2 23 ( 2 23 ⋅ E y + M y ) + μ − 127 ] \frac1{2^{23}}{\left(\textcolor{green}{2^{23}\cdot E_S+M_S}\right)}+\mu-127=-\frac12\left[\frac1{2^{23}}{\left(\textcolor{red}{2^{23}\cdot E_y+M_y}\right)}+\mu-127\right] 2 2 3 1 ( 2 2 3 ⋅ E S + M S ) + μ − 1 2 7 = − 2 1 [ 2 2 3 1 ( 2 2 3 ⋅ E y + M y ) + μ − 1 2 7 ]
化简,得到
2 23 ⋅ E S + M S = 3 2 ⋅ 2 23 ( 127 − μ ) − 1 2 ( 2 23 ⋅ E y + M y ) \textcolor{green}{2^{23}\cdot E_S+M_S}=\frac32\cdot 2^{23}(127-\mu)-\frac12 \left(\textcolor{red}{2^{23}\cdot E_y+M_y}\right) 2 2 3 ⋅ E S + M S = 2 3 ⋅ 2 2 3 ( 1 2 7 − μ ) − 2 1 ( 2 2 3 ⋅ E y + M y )
上面说过y = 2 23 ⋅ E y + M y , S = 2 23 ⋅ E S + M S y=\textcolor{red}{2^{23}\cdot E_y+M_y},S=\textcolor{green}{2^{23}\cdot E_S+M_S} y = 2 2 3 ⋅ E y + M y , S = 2 2 3 ⋅ E S + M S
S = 3 2 ⋅ 2 23 ( 127 − μ ) − y 2 \textcolor{green}{S}=\frac32\cdot 2^{23}(127-\mu)- \textcolor{red}{\frac y2} S = 2 3 ⋅ 2 2 3 ( 1 2 7 − μ ) − 2 y
所以S = n u m b e r − y 2 S=number-\frac y2 S = n u m b e r − 2 y n u m b e r number n u m b e r
实际上该常数的由来目前还是未知的,David Eberly博士曾经做出了最接近的解释,但实际上是错误的.不同的论文给出了不同的取值,并且都有各自的理由,而且一致认为源代码里的值是不好的,这里不做讨论,直接使用源代码的值代入
1 i = 0x5F3759DF - (i >> 1 );
由于float无法进行位运算,所以将它转换成long后再位移,虽然这样会损失精度,但是影响可以忽略不计.
此时已经运算完成,再把long转换回float
当前得到的y y y
设y y y x x x
1 + x = 1 ln 2 1+x=\frac{1}{\ln 2} 1 + x = ln 2 1
转换为x关于y的函数,得到
y = 1 x y=\frac{1}{\sqrt x} y = x 1
利用牛顿迭代法
f ( y ) = 1 y 2 − x = 0 f(y)=\frac 1{y^2}-x=0 f ( y ) = y 2 1 − x = 0
带入X n = y X_n=y X n = y
x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1}=x_n -\frac{f(x_n)}{f^{'}(x_n)} x n + 1 = x n − f ′ ( x n ) f ( x n )
化简
y n + 1 = y − f ( y ) f ′ ( y ) = y − 1 y 2 − x − 2 y 3 = y + y − y 3 ⋅ x 2 \begin{array}{ll} y_{n+1} &=y-\frac{f(y)}{f^{'}(y)} \\\\ &=y-\frac{ \frac{1}{y^2}-x }{ -\frac{2}{y^3} } \\\\ &=y+\frac{y-y^3 \cdot x}{2} \end{array} y n + 1 = y − f ′ ( y ) f ( y ) = y − − y 3 2 y 2 1 − x = y + 2 y − y 3 ⋅ x
得到最后一行代码.
1 y = y * (threehalfs - (x2 * y * y));
