1. 《雷神之锤 III 竞技场》快速求平方根导数源代码
- 此处先列出其源代码,这段代码的目标是计算一个浮点数平方根的导数,也就是如下形式: f ( x ) = 1 x f(x) = \frac{1}{\sqrt{x}} f(x)=x1
- 这段代码可以说非常难以理解,尤其是 0 x 5 f 3759 d f 0x5f3759df 0x5f3759df 更是莫名奇妙,不知所云。
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 // 2nd iteration, this can be removed// y = y * ( threehalfs - ( x2 * y * y ) ); return y;
}
2. 前置知识
- IEEE754-浮点数的表示方法,可以参考我的另外一篇文章 《IEEE754 -浮点数的表示》。
3. 推导过程
3.1 第一步:对要计算的公式两边同时取 l o g log log
-
正本清源,我们要求的就是这样一公式 y = 1 x y = x − 1 2 y = \frac{1}{\sqrt{x}} \\ y = x^{- \frac{1}{2}} y=x1y=x−21
-
对等式两边同时取以2为底的对数
不要问为什么取对数,取对数并不是目的,这只是一种 计算的技巧,在这里只需要知道等式两部同时取对数后仍然是相等的
log 2 ( y ) = log 2 ( x − 1 2 ) log 2 ( y ) = − 1 2 log 2 ( x ) \log_2(y) = \log_2(x^{-\frac{1}{2}}) \\ \log_2(y) = -\frac12\log_2(x) log2(y)=log2(x−21)log2(y)=−21log2(x)
3.2 第二部:计算出对浮点数取log的整数表达形式
-
将《IEEE754 -浮点数的表示》 中提到的浮点数表达式带入 第一步 所推导出的公式, 计算出对浮点数取 l o g log log的整数表达形式
-
《IEEE754 -浮点数的表示》 提到浮点数的计算公式为: ( − 1 ) S ∗ ( 1 + F 2 23 ) ∗ 2 ( E − 127 ) (-1)^S*(1+\frac F{2^{23}})*2^{(E-127)} (−1)S∗(1+223F)∗2(E−127)
-
不考虑符号位。对其取对数
log 2 ( ( 1 + F 2 23 ) ∗ 2 E − 127 ) \log_2\left(\left(1+\frac{\mathrm{F}}{2^{23}}\right)*2^{\mathrm{E}-127}\right) log2((1+223F)∗2E−127) -
化简一步
log 2 ( 1 + F 2 23 ) + E − 127 \log_2\left(1+\frac{\mathrm{F}}{2^{23}}\right)\:+\:\mathrm{E}\mathrm{-127} log2(1+223F)+E−127 -
推导到这一步似乎没办法化简了,我们知道 F 2 23 \frac{F}{2^{23}} 223F 是一个大于0,小于1的数,令 x = F 2 23 x= \frac{F}{2^{23}} x=223F,我们画出 log 2 ( 1 + x ) \log_2(1+x) log2(1+x)的曲线图,在0 - 1的范围内近乎是一条直线,所以 log 2 ( 1 + x ) ≈ x \log_2(1+x) \approx x log2(1+x)≈x
-
加入误差修正(假设误差为 μ μ μ),如下图修正0.045的效果看起来更好一点,所以 log 2 ( 1 + x ) ≈ x + μ \log_2(1+x) \approx x+μ log2(1+x)≈x+μ
此处的关键就是找到一个 μ μ μ, 使得总体误差最小,根据前人的计算 μ = 0.0450465 μ = 0.0450465 μ=0.0450465左右,整体误差最小,
-
带入上面公式
log 2 ( 1 + F 2 23 ) + E − 127 ≈ x + μ + E − 127 ≈ F 2 23 + μ + E − 127 ≈ F 2 23 + E + u − 127 ≈ F 2 23 + 2 23 E 2 23 + u − 127 − − − − − − ≈ 1 2 23 ( F + 2 23 ∗ E ) + u − 127 \log_2\left(1+\frac{\mathrm{F}}{2^{23}}\right)\:+\:\mathrm{E}\mathrm{-127}\\ \approx x+μ+E -127 \\ \approx \frac{F}{2^{23}}+μ+E -127 \\ \approx \frac{F}{2^{23}}+E+u -127 \\ \approx \frac{F}{2^{23}}+\frac{2^{23}E}{2^{23}}+u -127 \\ ------\\ \approx \frac{1}{2^{23}} (F +2^{23}*E)+u -127 log2(1+223F)+E−127≈x+μ+E−127≈223F+μ+E−127≈223F+E+u−127≈223F+223223E+u−127−−−−−−≈2231(F+223∗E)+u−127 -
观察最终计算出的公式,红框中的部分就是 《IEEE754 -浮点数的表示》所提到的 浮点数二进制的整形表达形式
3.3 第三部:将 第二步 带入 第一步
- 化简一下公式
l o g 2 ( y ) = − 1 2 log 2 ( x ) 1 2 23 ( F y + 2 23 ∗ E y ) + μ − 127 = − 1 2 ( 1 2 23 ( F x + 2 23 ∗ E x ) + μ − 127 ) − − − − − − ( F y + 2 23 ∗ E y ) = 3 2 ∗ 2 23 ( 127 − μ ) − 1 2 ( F x + 2 23 ∗ E x ) log_2(y) = -\frac12\log_2(x) \\ \frac{1}{2^{23}}(\mathrm{F}_{y}+2^{23}*\mathrm{E}_{y})+\mu-127=-\frac{1}{2}\left(\frac{1}{2^{23}}(\mathrm{F}_{x}+2^{23}*\mathrm{E}_{x})+\mu-127\right) \\ ------\\ (\mathrm{F}_y+2^{23}*\mathrm{E}_y)=\frac32*2^{23}(127-\mu)-\frac12(\mathrm{F}_x+2^{23}*\mathrm{E}_x) log2(y)=−21log2(x)2231(Fy+223∗Ey)+μ−127=−21(2231(Fx+223∗Ex)+μ−127)−−−−−−(Fy+223∗Ey)=23∗223(127−μ)−21(Fx+223∗Ex) - 结果已经呼之欲出了,如下图所示,红框代表的是 y y y的二进制值, 蓝框代表的是 x x x的二进制值。
y 二进制 = 3 2 ∗ 2 23 ( 127 − μ ) − 1 2 ( x 二进制 ) 令 μ = 0.0450465 (这其实是根据 0 x 5 f 3759 d f 反推出来的),使误差最小 − − − − y 二进制 = 0 x 5 f 3759 d f − 1 2 ( x 二进制 ) y_{二进制} = \frac32*2^{23}(127-\mu)-\frac12(x_{二进制})\\ 令μ = 0.0450465(这其实是根据0x5f3759df反推出来的),使误差最小 \\ ----\\ y_{二进制} = 0x5f3759df-\frac12(x_{二进制}) y二进制=23∗223(127−μ)−21(x二进制)令μ=0.0450465(这其实是根据0x5f3759df反推出来的),使误差最小−−−−y二进制=0x5f3759df−21(x二进制)
4.最后在使用牛顿迭代法再次逼近一次结果,提高精度
- 经过上面推导,已经求出来一个值了,这个值存在一定的误差,如果接受不了这个误差,就可以使用牛顿求根法进行一次或多次逼近以提高精度,如果能够接受这种误差便可以结束了
4.1 认识牛顿迭代
-
牛顿迭代法的核心如下图所示,通俗点来讲就是通过已知点,不断逼近 f ( x ) = 0 f(x) = 0 f(x)=0的点
-
通用公式如下,很好理解: x n + 1 = x n − f ( x n ) f ’ ( x n ) x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)} xn+1=xn−f’(xn)f(xn)
以 x 0 x_0 x0到 x 1 x_1 x1为例,通过 y 1 y_1 y1除以斜率计算出 x 0 − x 1 x_0 - x_1 x0−x1 之间的距离,斜率就是在 y 0 y_0 y0 这个点上求导
4.2 对 y = 1 x y = \frac{1}{\sqrt{x}} y=x1进行牛顿迭代
- 两边同时求平方进行一步变换,目的是消除根号
y = 1 x y 2 = 1 x y = \frac{1}{\sqrt{x}} \\ y^2 = \frac{1}{x} y=x1y2=x1 - 两边同时取倒数,目的是将方程转换为以 y y y 为变量的表达式 1 y 2 = x \frac{1}{y^2} = x y21=x
- 整理一下,得到 f ( x ) = 0 f(x) = 0 f(x)=0 的形式: 1 y 2 − x = 0 \frac{1}{y^2} - x = 0 y21−x=0
- 回想一下我们的目标, x x x 是我们要开计算平方根导数的值, y y y是要计算的结果,我们之前已经计算出了一个不太精确的 y y y,现在的目标就是通过这个不太精确的 y y y逼近下一个点,
- 带入牛顿迭代法的通用公式 f ( y ) = 1 y 2 − x f ′ ( y ) = − 2 y 3 − − − − − − − − − − − y n + 1 = y n − f ( y n ) f ′ ( y n ) − − − − − − − − − − − y n + 1 = y n − 1 y n 2 − x − 2 y n 3 y n + 1 = y n ( 3 − x y n 2 ) 2 − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − y n + 1 = y n ( 1.5 − 0.5 x y n 2 ) f(y) = \frac{1}{y^2} - x\\f'(y) = \frac{-2}{y^3}\\ -----------\\ y_{n+1}=y_n-\frac{f(y_n)}{f'(y_n)}\\ -----------\\ y_{n+1}=y_n-\frac{\frac{1}{y_n^2}-x}{\frac{-2}{y_n^3}}\\ y_{n+1}=\frac{y_n(3-xy_n^2)}2\\ ---------------------------------\\ y_{n+1}=y_n(1.5-0.5xy_n^2) f(y)=y21−xf′(y)=y3−2−−−−−−−−−−−yn+1=yn−f′(yn)f(yn)−−−−−−−−−−−yn+1=yn−yn3−2yn21−xyn+1=2yn(3−xyn2)−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−yn+1=yn(1.5−0.5xyn2)
- 至此,公式与代码都对上了