首页 资讯 倒数速算法(FastRootRoot)[3],牛顿迭代法

倒数速算法(FastRootRoot)[3],牛顿迭代法

更新时间:2022-09-25 11:07:16 分类:资讯 浏览:25

简介

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

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

版权声明: 本站内容部分来源网络,版权归作者所有,如有侵权,请联系我们删除!