《微微一笑很倾城》中肖奈大神说的平方根倒数速算法是什么鬼？三十分钟理解！

最新推荐文章于 2024-05-12 11:23:43 发布

大饼博士X

最新推荐文章于 2024-05-12 11:23:43 发布

阅读量5.6k

点赞数 7

分类专栏：三十分钟理解系列文章标签：算法加速

本文链接：https://blog.csdn.net/xbinworld/article/details/52421847

版权

三十分钟理解系列专栏收录该内容

9 篇文章 17 订阅

订阅专栏

在电影《微微一笑很倾城》中，肖奈大神在玻璃上写了一堆公式，提到平方根倒数速算算法，这个到底是一个什么算法？笔者看电影的时候打开手机学了一下，发现该算法的作者真乃神人！今天有空，就把该算法写一写。

在3D图形编程中，经常要求平方根的倒数，即1/Sqrt(x)，如果用一般的代码(float)(1.0/sqrt(x)),,精度高，但是非常慢；我们需要一个快速，而又足够高精度的算法；著名游戏《雷神之锤III》的代码在2002年左右被披露，人们发现了一段用于快速计算平方根倒数的代码，下面是整理后的代码（去掉了一些宏定义）。

float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f3759df - (i >> 1);// L1：gives initial guess y0
x = *(float*)&i; //l2：convert bits back to float
x = x*(1.5f - xhalf*x*x);  //l3：Newton step, repeating increases accuracy
return x;
}

该程序运行效率极高，经测试，基本是使用直接开根号求倒数程序的4倍速度！一时间惊为天人。那么这段代码到底怎么理解？为什么中间出现了0x5f3759df这样一个完全无厘头的magic number？

--------------------------------------------------------------------------------------------------------------------------------------------------

该算法的本质其实就是牛顿迭代法（Newton-Raphson Method，简称NR）。NR是一种求方程的近似根的方法。首先要估计一个与方程的根比较靠近的数值，然后根据公式推算下一个更加近似的数值，不断重复直到可以获得满意的精度。其公式如下：

函数：y=f(x)
其一阶导数为：y'=f'(x)
则方程：f(x)=0 的第n+1个近似根为
x[n+1] = x[n] - f(x[n]) / f'(x[n])

--------------------------------------------------------------------------------------------------------------------------------------------------

求一个数a的平方根的倒数，实际就是求方程f(x)=1/(x^2)-a=0的解；将该方程按牛顿迭代法的公式展开为：

x[n+1]=x[n]*(3/2-a/2*x[n]*x[n])

也就是上面代码蓝色部分的第L3行；所以很明显，这段代码就是进行一次迭代的牛顿迭代法。只进行一次迭代就结束程序，又要保证精度，也就是说，牛顿迭代法的初始值要非常精准，才能在一次迭代后完成计算。整个代码最精彩的就是L1行了，i = 0x5f3759df - (i >> 1);到底是什么意思？为什么初始值可以这样算呢？了解该代码，需要先了解一些float point和fix point的表示方法[2]：

一个浮点数float是由32位二进制位表示的有理数，分为三部分。其中符号占1位，表示正负，记为Si；指数占接下来的8位，表示经过偏移处理后的指数，即实际表示E（如图中为124），需要偏移B（图中为2的8次方减1,127。B为一个固定值），最后得指数值为E-B；有效数字（除最高位1以外，因为前面有指数，所以一个数肯定可以表示成1.xxxxx * 2^(kkk)，只保留xxxxx）占剩下的23位，记为m（0<m<1）,图中的 $\scriptstyle m=1\times 2^{-2}=0.250$ 。

所以浮点数的结构公式为： $\scriptstyle x=(-1)^{\mathrm{Si}}\cdot(1+m)\cdot 2^{(E-B)}$ ，图中 $\scriptstyle x=(1+0.250)\cdot 2^{-3}=0.15625$

整数fix point的表示相对简单，符号占1位，数值占剩下的31位。如果用上图的浮点数字节序列来表示整数，那么 $\scriptstyle I=E\times 2^{23}+M$ ，即 $I=124\times 2^{23} + 2^{21}$ ；平方根倒数函数仅能处理正数，所以符号位均为0。

对于同样的32位二进制数码，若为浮点数表示时实际数值为 $\scriptstyle x=(1+m_x)2^{e_x}$ ，而若为整数表示时实际数值则为 $\scriptstyle I_x=E_xL+M_x$ ，其中 $\scriptstyle L=2^{n-1-b}$ ，这里n=32，b=8。式子中引入的新变量为：

$m_x=\frac{M_x}{L}$ ------------------------------------等式1

$e_x=E_x-B$ ，其中 $B=2^{b-1}-1$ ------------等式2

理解浮点数和整数的表示后，下面开始推导。

x的平方根倒数方程为：

$y=\frac{1}{\sqrt{x}}$

两边取对数有：

$\log_2{(y)}=-\frac{1}{2}\log_2{(x)}$

因为浮点数可表示为： $\scriptstyle x=(1+m_x)2^{e_x}$ ，所以也有，代入上式有：

$\log_2(1+m_y)+e_y=-\frac{1}{2}\log_2{(1+m_x)}-\frac{1}{2}e_x$

由于 $\scriptstyle 0 \le x < 1$ ，有 $\scriptstyle \log_2{(1+x)}\approx {x}$ ，则在此可定义 $\sigma$ 与x的关系为 $\scriptstyle \log_2{(1+x)}\cong x+\sigma$ ， $\sigma$ 表示误差，所以将 $\scriptstyle \log_2{(1+x)}= x+\sigma$ 代入上式得：

$m_y+\sigma+e_y=-\frac{1}{2}m_x-\frac{1}{2}\sigma-\frac{1}{2}e_x$

将等式1，等式2代入到上述方程中，有：

$M_y+(E_y-B)L=-\frac{3}{2}\sigma{L}-\frac{1}{2}M_x-\frac{1}{2}(E_x-B)L$

移项整理得：

$E_yL+M_y=\frac{3}{2}(B-\sigma)L-\frac{1}{2}(E_xL+M_x)$

又因为浮点规格存储的正浮点数x，若将其作为整数表示，则示值为： $\scriptstyle I_x=E_xL+M_x$ ，所以x的平方根倒数的首次近似值的整数表示值为：

$I_y=E_yL+M_y=R-\frac{1}{2}(E_xL+M_x)=R-\frac{1}{2}I_x$ ，

其中 $R=\frac{3}{2}(B-\sigma)L$ ， $B=2^{b-1}-1$ ， $\scriptstyle L=2^{n-1-b}$ ，n=32，b=8。

这个式子对应着源代码中的这一行：i = 0x5f3759df - ( i >> 1 );，然后将整数表示值换回表示浮点数：x = * ( float * ) &i;。这样就得到了浮点数的平方根倒数的近似值。

0x5f3759df 对应着R，即3/2(B- $\sigma$ )L.当R为0x5f3759df时，有 $\scriptstyle \sigma=0.0450461875791687011756$ .

"现在不仅该算法的原作者不明，人们也仍无法明确当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验：他编写了一个函数，以在一个范围内遍历选取R值的方式将逼近误差降到最小，以此方法他计算出了线性近似的最优R值0x5f37642f（与代码中使用的0x5f3759df相当接近），但以之代入算法计算并进行一次牛顿迭代后，所得近似值与代入0x5f3759df的结果相比精度却仍略微更低。……在Charles McEniry的论文中，他使用了一种类似Lomont但更复杂的方法来优化R值：他最开始使用穷举搜索，所得结果与Lomont相同；而后他尝试用带权二分法寻找最优值，所得结果恰是代码中所使用的魔术数字0x5f3759df"---维基百科

参考资料：

[1] CHRIS LOMONT, FAST INVERSE SQUARE ROOT

[2] http://blog.csdn.net/heloowird/article/details/21862251

[3] http://blog.chinaunix.net/uid-9255716-id-107951.html

[4] http://blog.csdn.net/xiaoguohaha/article/details/21652643