倒数速算法(FastRootRoot)[3],牛顿迭代法
简介
Fast Root[3] 是一个相当有名的算法,它的具体作者不再可用。因其在《雷神之锤》中广泛使用该算法而闻名。
为什么需要平方根的倒数而不是平方根本身?因为图形中大量的运算需要用到平方根的倒数,而不是平方根,比如frac{1}{sqrt{x^2 + y^2 +的()运算z^2}}(x,y, z) 包括倒数平方根;并且图形中每一帧的渲染至少需要数万次调用。如果能加快倒数平方根的计算,渲染速度会大大提高。
算法诞生于 20 多年前。今天,该算法只有研究价值,没有使用价值[1]。 CPU内置指令,速度更快平方根公式,错误更少;我们只需要使用 1/sqrt (y)。
不过,这并不影响我们对这个算法的精妙之处的理解。仅仅4行C代码,就涉及到1)IEEE 754标准2)C语言3)牛顿迭代法。
这里的重点是第二步,也就是如何衍生魔法。详情请参考[2] [4]。
算法
//from wikipedia [3]
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
/*核心代码*/
// 1. reinterpret_cast from float to int
i = * ( long * ) &y; // evil floating point bit level hacking(对浮点数的邪恶位元hack)
// 2. 估算平方根倒数
i = 0x5f3759df - ( i >> 1 ); // what the fuck? (....why not 0x69696969?)
// reinterpret_cast from int to float
y = * ( float * ) &i;
// 3. 牛顿法
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
算法大致分为三部分:
求平方根导数
在 IEEE 754 标准中,一个 32 位浮点数包含一个符号位、一个 8 位指数 E 和一个 23 位尾数 M。任何浮点数都可以表示为:(1.M) times 2^{ E - 127} = (1 + frac{M}{2^{23}}) times 2^{E - 127}, M in[0, 2^{ 23} ), E in [0, 256) .
目的:需要计算yfrac{1}{sqrt{y}} = y^{-frac{1}{2}}的倒数平方根,y的IEEE 754形式为y = (1 + frac{M_y}{2^{23}}) * 2^{E_y - 127}。
为了数学直觉,我们取 y 的以 2 为底的对数:
begin{align}log y &= log{((1 + M_y) * 2^{E_y - 127})} \&= log{(1+frac{M_y}{2^ {23}})} + E_y - 127end{align}
frac{M_y}{2^{23}}显然在区间[0,1]内,而在这个区间内,有log(1+x)x;
图片-
如果加上常数 mu,近似值会更好:
图片-
使用这个近似值,可以进一步简化:
begin{align}log y&= log{(1+frac{M_y}{2^{23}})} + E_y - 127 \& frac{M_y}{2^{23 }} + mu + E_y - 127 \&= frac{1}{2^{23}}(M_y + E_y * 2^{23}) + mu - 127end{align}
而M_y + E_y * 2^{23}正是float类型y转换为int类型后的值。
现在对 y^{-frac{1}{2}} 取相同的对数:
log y^{-frac{1}{2}} = -frac{1}{2} log y -frac{1}{2}left( frac{1}{ 2^{23}}(M_y + E_y * 2^{23}) + mu - 127 right)
假设y^{-frac{1}{2}}的准确值为A,其指数部分和尾数部分分别为E_A和M_A,同样取A的对数:
log A frac{1}{2^{23}}(M_A + E_A * 2^{23}) + mu - 127
上面两个方程用不同的方式表示$log A $,即log y^{-frac{1}{2}},所以上面两个方程可以解得到M_A + E_A * 2^ {23}:
begin{align}log A &= log y^{-frac{1}{2}}\-frac{1}{2}left( frac{1}{2^ {23}}(M_y + E_y * 2^{23}) + mu - 127 right) &= frac{1}{2^{23}}(M_A + E_A * 2^{23}) + 亩 - 127 \end{对齐}
省略中间过程得到:
begin{align}M_A + E_A * 2^{23} &= 2^{23} * frac{3}{2}(127 - mu) - frac{1}{2}(M_y + E_y * 2^{23}) \end{align}
当第一项中的mu取$0.$[3]时,第一项的结果是,十六进制表示是;如果第二项的frac{1}{2}变成了移位操作,上面的公式和代码中的这一行很相似
i = 0x5f3759df - ( i >> 1 );
这有点离奇平方根公式,如果 mu 的值是正确的,你应该得到完全相同的数字。
牛顿迭代法
牛顿迭代法的更新公式为:
x_{t+1} = x_t - f(x_t) / f^prime(x_t)
牛顿迭代法的前提条件:要求初始点在函数的零点附近,因为 i = - ( i >> 1 );估计结果误差应该很小,可以认为前提条件得到满足。
此时,我们的问题是,给定 y,我们需要找到 x 使得 frac{1}{sqrt{y}} = x,即 frac{1}{x^2} = y,移项 Get f(x) = frac{1}{x^2} - y = 0 ,对 f(x) 应用牛顿迭代得到迭代公式:
begin{align}x_{t+1} &= x_t - f(x_t) / f^prime(x_t) \x_{t+1} &= frac{3}{2}x_t - frac{y}{2}x_t^3end{align}
翻译成代码是:
y = y * ( threehalfs - ( x2 * y * y ) );
参考文献
[1] 还在使用快速 root 吗? /Is-fast---root-still-used
[2] Fast Root — A Quake III /watch?v=
[3] wiki-倒数平方根算法/zh-hans/%E5%B9%B3%E6%96%B9%E6%A0%B9%E5%80%92%E6%95%B0%E9% 80%9F%E7%AE%97%E6%B3%95
[4]论文//2003/.pdf
本文是使用知乎上创建和发布的