快速平方根取倒数算法

很久以前就想写一下这个快速平方根取倒数算法了，这个好像还很有历史“渊源”来着，据说给当时“雷神之锤”的开发带来了巨大优化。

首先我们思考一下，平方根应该用哪种算法求，显然这是一个只能逼近的值，因此我们想到了“牛顿迭代”算法：

\[x_n=x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})} \]

考虑函数 \(f(x)=\frac{1}{x^2}-a\) ，先选定一个初始值 \(x_0\) 然后：

\[\begin{aligned}x_n&=x_{n-1}-\frac{ax_{n-1}^3-x_{n-1}}{2}\\ &=\frac{3x_{n-1}-ax_{n-1}^3}{2}\end{aligned} \]

迭代求解。

显然复杂度是与求解的精度相关的，迭代次数和答案的精度正相关。那么，能否选定某个恰当的 \(x_0\)，使得迭代次数得到大幅的减少呢？答案是肯定的。

早期游戏工程师为了简化复杂度（众所周知早期算法工程师常从底层最大限度的优化算法，极致压榨 \(PC\) 机效率），从浮点数的性质入手：

\(32\) 位浮点数的 IEEE754 存储形式是这样的:

• \(S\) 是符号位，在开平方根是不必理会，因为计算机系统运算是基于实数的，在实数中对负数开平方是没有意义的
• \(E\) 表示阶码，但是有 \(127\) 的平移用来表示正负，故 \(E−127\) 是阶码表示的实际值
• \(M\) 表示尾码，表示原二进制浮点数第一个 \(1\) 之后的所有数位，基本可以认为隐藏了一位 \(1\) ,实际值为 \(1.xxxxxxx\)。

那么类似于十进制的科学计数法，我们可以将浮点数表示为：

\[y=(1+ \frac{M}{2^{23}})⋅2^{E−127} \]

对于浮点数 \(x=\frac{1}{\sqrt{a}}\)，我们不妨取对数；

\[\begin{aligned} -\frac{1}{2}\log_2a&=-\frac{1}{2}[\log_2(1+\frac{M_a}{2^{23}})+E_a-127]\end{aligned} \]

我们观察到 \(\frac{M_a}{2^{23}}\in(0,1)\)，故 \(\log_a(1+x)\) 可近似为 \(x+\epsilon\) （\(\epsilon\) 为修正值），即：

\[\begin{aligned}\log_2\frac{1}{\sqrt{a}}&=-\frac{1}{2}(\frac{M_a}{2^{23}}+E_a+\epsilon-127)\\\log_2\frac{1}{\sqrt{a}}&=-\frac{1}{2}(\frac{M_a+2^{23}E_a}{2^{23}}+\epsilon-127)\\ \end{aligned} \]

观察这个部分 \(M_a+2^{23}E_a\)，非常巧妙地，跟据浮点数的存储规则，我们发现这个部分就是浮点数 \(a\) 强制转换为整数型后的表示，我们把这个部分称为 \(X_a\)。接下来，我们将 \(x\) 和 \(a\) 分别表示并联立：

\[\left\{ \begin{aligned}\log_2\frac{1}{\sqrt{a}}&=-\frac{1}{2}(\frac{X_a}{2^{23}}+\epsilon-127)\\ \log_2x&=\frac{X_x}{2^{23}}+\epsilon-127 \end{aligned} \right. \]

化简得；

\[\begin{aligned}\frac{X_x}{2^{23}}+\epsilon-127&=-\frac{1}{2}(\frac{X_a}{2^{23}}+\epsilon-127)\\X_x&=(381+\epsilon)2^{23}-\frac{1}{2}X_a \end{aligned} \]

接下来，为了使修正值更加正确，我们不妨提出一种启发式的求解方法：使 \(\log2(1+x)\) 与 \(x+\epsilon\) 在区间 \([0,1]\) 上平均差值为 \(0\)，即：

\[\begin{aligned}\int_0^1\log_2(1+x)dx-(x+\epsilon)dx&=0\\ \int_0^1\log_2(1+x)dx-\int_0^1(x+\epsilon)dx&=0\\ \int_0^1\log_2(1+x)dx&=\int_0^1(x+\epsilon)dx\\ \frac{1}{\ln2}\int_1^2\ln{u}du&=\frac{1}{2}+\epsilon\\ \frac{1}{\ln2}(2ln{2}-1)&=\frac{1}{2}+\epsilon\\ 2-\frac{1}{ln{2}}&=\frac{1}{2}+\epsilon\\ \epsilon&=\frac{3}{2}-\frac{1}{\ln2}\\\end{aligned} \]

最后将\(2^{22}\times(381-\epsilon)\) 十六进制硬编码为 0x5f3c551d，而 0x5f3759df 这个经典算法中的”魔数“怎么来的希腊奶（事实上确实是经典算法中的“魔数”能带来更加精确的答案。我猜是暴力二分试出来的，毕竟是程序员又不是数学家，费这劲干啥）。

综上，我们就完整的得到了精度非常高的初始值 \(x_0\)，接着将这种初始值代入牛顿法仅需 \(1\) 次迭代就能达到相当高的精度，大大降低了复杂度。

上个原版代码；

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 hackingi  = 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 removedreturn y;
}