倒数速算法(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

本文是使用知乎上创建和发布的