2021-01-23 快速平方根求倒数算法

此算法很有名，来自《雷神之锤3》Quake引擎源码q_math.c的552行，传言是3d游戏开发大神约翰卡马克写的，代码出来之后很震惊，不过据考证最早的实现人并不是他，至于是谁第一个实现的至今成迷，源码如下：

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 );               // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration

    return y;
}

这篇文章内容不是原创，主要来自网上一些介绍，博主一开始也没看明白这算法，看了一些网上介绍大概了解了，这里分析一下算法原理，同时做个备忘录：

算法有两步是核心，一个就是魔数：0x5f3759df,还有一个就是
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration

下面介绍如下：
输入 $a$ , 目标是计算 $y=\sqrt{1/a}$ ,列出方程就是:
$1/y^2 -a =0$
令 $f(y)= 1/y^2 - a$ 则目标就是求解这个方程的根,
这里有牛顿迭代法可以求解，牛顿迭代介绍比较多，这里不再缀述，直接给出公式:
$y_{n+1}=y_{n}-f(y_{n})/f'(y_{n})$
将 $y=1/y^2-a$ 代入上式并化简，得到:
$y_{n+1}=y_{n}*(1.5-0.5a*y_{n}*y_{n})$
这其实就是对应代码中的:

y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration

由此，可见代码中只进行了一次牛顿迭代就默认收敛了，那么为何作者那么确定能快速收敛，这主要取决于迭代初始值的选择，如果初始值与目标值越近，收敛越快

任何一个浮点数，在计算机中存储如下：
$s e⋯em⋯m$
s就1位，符号位置， e...e共8位 m...m共23位取值全部为0,1
如果用2进制来看（下文中默认都是2进制）
令 $0.m...m$ 的值为 $m'$ , $e' = e...e - B$
(浮点数为了包含负指数情形，因此上式有个偏移量 $B=127$ )
则上述浮点数的值为:
$(1+m')\times 2^{e'}$ 并且 $m'$ 取值范围在 $[0,1]$ 之间

如果把这个浮点数看做一个整数（也就是直接强转为整数,下同），并且令
$M =m...m$ $E=e...e$ $L=2^{23}$ 则这个整数取值为:
$M\times L+E$ （公式1）

同时 $M,m',E,e'$ 之间存在如下关系:
$m'=M/L L=2^{23} ,e'=e-B,B=127$

现在先把这些公式放到一边，我们回到原问题，其实要求一个迭代初始值，这个初始值是这么推导来的:
要计算 $y=1/sqrt(a)$ 两边取对数，得到:
$log_2{y}=-1/2 log_2{a}$ 现在用上面的浮点数表示替换, $y,a$ 得到
$log_2(1+m_y')+e_y' = -1/2[ log_2(1+m_a')+e_a']$
可以看见两边都有形如 $log_2(1+X)（0<=X<=1）$ 的形式
如果把 $log_2(1+X) (0<=X<=1)$ 这个函数近似为线性函数: $X+\sigma$ 代入上式，得到:
$m_y'+\sigma+e_y'=−1/2(m_a'+\sigma+e_a')$
再用上面的公式，可以简化为:
$3/2L(\sigma−B)+M_y+LE_y=−1/2(M_a+LE_a)$
如果记:
$I_y=M_y\times L+E_y$
$I_a =M_a\times L +E_a$
$k =3/2L(B-\sigma)$
有:
$I_y=(-1/2)I_a + k$ (2)
其中 $k$ 对应的就是魔数: 0x5f3759df ，

根据上文分析，表达式(2)意味着：
根据公式(1), 如果把输入的浮点数a,看做整数，其取值就是 $I_a$
则根据(2)计算出的就是浮点数 $y=1/sqrt(a)$ 看成整数的值，
再把这个值强转为浮点数，就是浮点数y的值，
(当然在计算过程中我们近似利用了 $log_2(X+1) = X+\sigma$
所以这样的算出的值是有误差的，但是误差比较小，作为牛顿迭代的初始值，只需要一次迭代就收敛了；)

而公式2对应的源码就是(已添加注释）:

    y  = number; //number相当于公式中的a
    i  = * ( long * ) &y;//a强转为整数，得到I_a             
    i  = 0x5f3759df - ( i >> 1 );    //i>>1就等于I_a/2 ,
    y  = * ( float * ) &i; //最后再强转为整数

据此，就分析完了，从上面可以看出魔数主要由 $\sigma$ 确定，至于 $\sigma$ 怎么确定的，感兴趣的读者可以深究，但是 $\sigma$ 的确定只与1个确定函数有关，即: $log_2(1+X)$ ,它的取值与算法的输入输出是无关的，也就是说，不管怎么取，只会影响魔数:0x5f3759df；而不会影响代码的结构

因此，对于我们的目的已经达到了，这个算法的原理已经解释清楚了

2021-01-23 快速平方根求倒数算法

2021-01-23 快速平方根求倒数算法

相关阅读更多精彩内容

友情链接更多精彩内容