计算机与数学 —— 雷神之锤3源码中的快速逆平方根算法

这篇博客介绍了在雷神之锤3源代码中快速求逆平方根的算法。


源码

雷神之锤3中的逆平方根算法如下:

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

函数接受一个浮点数作为输入,输出的结果是平方根的倒数。

单精度浮点数

在计算机中,单精度浮点数使用32位来储存表示:

Floating Point

在这32位中,最高位为符号位,后面的8位为整数ex,代表浮点数的指数,而最后的23位表示的是小数部分mx,小数部分第一位表示 2 − 1 2^{-1} 21,第二位表示 2 − 2 2^{-2} 22……

所以如果x是一个正浮点数,则有:

x = 2 e x ( 1 + m x ) x = 2^{e_{x}} (1+m_{x}) x=2ex(1+mx)

另外,如果需要把一个浮点数转化为整数形式,则需要做如下的运算:

I x = E X L + M x = L ( e x + B + m x ) Ix = E_XL + M_x=L(e_x+B+m_x) Ix=EXL+Mx=L(ex+B+mx)

在上面的公式中,L指的是指数部分需要的唯一次数,这里是 2 23 2^{23} 223,B是127,而M是小数部分对应的整数版本。

牛顿法

牛顿迭代法(简称牛顿法),是用于计算机求解任意连续函数的根值的一种方法。

假设对于如下的函数,我们想要求这个函数的根:

Function

该如何求得这个根呢?首先,我们先猜一个解 x 0 x_{0} x0,并且认为它是函数的解。

但是由于它其实并不是函数的解,那么现在,我们需要将这个解进行迭代,从而让其逼近真正的解。

在此处,我们可以在 ( x 0 , y 0 ) (x_{0}, y_{0}) (x0,y0)处作其切线,求得该直线的方程:

y − f ( x 0 ) = f ′ ( x 0 ) ( x − x 0 ) y - f(x_0) = f'(x_0)(x - x_0) yf(x0)=f(x0)(xx0)

并且求直线的根,此时会发现已经对于真正的解逼近了一步。

推广到n,继续迭代,就可以足够逼近真正的解了:

x i + 1 = x i − f ( x i ) f ′ ( x i ) x_{i+1} = x_{i} - \dfrac{f(x_i)}{f'(x_i)} xi+1=xif(xi)f(xi)

此时发现 f ( x i ) f(x_i) f(xi) f ′ ( x i ) f'(x_i) f(xi)可以被一个统一的函数 g ( x i ) g(x_i) g(xi)来表示:

g ( x ) = f ( x ) f ′ ( x ) g(x) = \dfrac{f(x)}{f'(x)} g(x)=f(x)f(x)

令ε为当前的解与真正解r的距离:
ϵ i = x i − r \epsilon_i = x_i - r ϵi=xir

综合上面三个方程,可得:

ϵ i + 1 = ϵ i − g ( x i ) \epsilon_{i+1} = \epsilon_i - g(x_i) ϵi+1=ϵig(xi)

因此只要ε值小于某个特定的值,我们可以认为此时的x和方程的解已经很接近了。

算法分析

如果需要求得一个浮点数的平方根倒数,方程如下:

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

转化为关于y的方程,有:

f ( y ) = 1 y 2 − x = 0 f(y) = \dfrac{1}{y^2} - x = 0 f(y)=y21x=0

转化为牛顿法使用的方程,有:

y n + 1 = y n ( 3 − x y n 2 ) 2 y_{n+1} = \dfrac{y_n(3 - xy_n^2)}{2} yn+1=2yn(3xyn2)

此时,对原本的方程等号两边同时取2的对数,就有:

log ⁡ 2 y = − 1 2 log ⁡ 2 ( 1 + m x ) \log_{2}{y} = -\dfrac{1}{2}\log_{2}{(1+m_x)} log2y=21log2(1+mx)

因为 m x ≥ 0 m_x\ge0 mx0并且 m x < 1 m_x<1 mx<1,那么在这个区间内,可以取近似为:
log ⁡ 2 ( 1 + m x ) ≈ m x + σ \log_2(1+m_x)\approx m_x + \sigma log2(1+mx)mx+σ

根据方差的计算,当 σ = 0.0430357 \sigma=0.0430357 σ=0.0430357时,整体的偏差是最小的,此时上面的等号两边应该相当。

因此,把上面的完全整合起来,最终的 I x I_x Ix可以写成:
L log ⁡ 2 x + L ( B − σ ) L\log_2{x}+L(B-\sigma) Llog2x+L(Bσ)

I y ≈ 3 2 L ( B − σ ) − 1 2 I x I_y\approx \dfrac{3}{2}L(B-\sigma)-\dfrac{1}{2}I_x Iy23L(Bσ)21Ix

最终,写成代码就是:

i = 0x5f3759df - ( i >> 1 );

3 2 L ( B − σ ) = 0 x 5 f 3759 d f \dfrac{3}{2}L(B-\sigma) = 0x5f3759df 23L(Bσ)=0x5f3759df

<全文完>

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页