Q_rsqrt代码阅读

代码


xxxxxxxxxx
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 hacking 
    i = 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 removed 
    return y; 
}

背景

第一眼看代码我们可以get到这个函数需要输入一个浮点数，输出一个浮点数，那我们接下来先拆解下这个输入输出的秘密。首先，我们知道一个单位向量的公式。

\begin{align*} \parallel v \parallel &= \sqrt{x^2 + y^2 + z^2} \\ \hat v &= \frac {v} {\parallel v \parallel} \\ \end{align*}

出于计算机硬件设计的历史原因，硬件除法指令和开方的实现性能较差，你可以查看intel的指令手册关于latency的数据。那我们可以从软件层规避掉这种指令吗？答案是可以的，我们接下来看看id software的工程师是如何处理这个问题的。

\begin{align*} \hat v &= v \times \frac {1} {\parallel v \parallel} \\ \frac {1} {\parallel v \parallel} &= \frac {1} {\sqrt{x^2 + y^2 + z^2}} \end{align*}

$\frac {1} {\parallel v \parallel}$ 。

y = \frac {1} {\sqrt x}

$x^2 + y^2 + z^2$ $y = \frac {1} {\sqrt {x^2 + y^2 + z^2}}$ 。

注释1: evil floating point bit level hacking

因为我们不能在浮点数上进行位运算的操作，我们先将float转成int。


xxxxxxxxxx
i = * ( long * ) &y; // evil floating point bit level hacking

可以看看下面quake 2 engine的一个函数Q_fabs(float absolute)。


xxxxxxxxxx
// from quake 2 code
// *(TYPE*) &VAR, since we can't perform bitwise operation on float point number
float Q_fabs (float f)
{
    int tmp = * ( int * ) &f;
    // set the sign bit
    // 7F => 0111 1111
    tmp &= 0x7FFFFFFF;
    return * ( float * ) &tmp;
}

注释2: what the fuck


xxxxxxxxxx
i = 0x5f3759df - ( i >> 1 ); // what the fuck?

$log_2$ 函数等式消除除法指令。

\begin{align*} y &= \frac {1} {\sqrt x} \\ log_2(y) &= log_2({x ^ {-\frac{1}{2}}}) \\ log_2(y) &= -\frac {1}{2}log_2(x) \tag{1} \end{align*}

$log_2(x)$ $3.75_{10}$ 十进制转二进制计算器 $11.11_{2}$ 。

11.11 = 1 \times 2 ^ {1} + 1 \times 2 ^ {0} + 1 \times 2 ^ {-1} + 1 \times 2 ^ {-2}

$d_{n-1}×r ^ {n-1} + d_{n-2}×r ^ {n-2} + ... + d_1×r^1 + d_0×r^0$ 计算十进制数值。

1.111 \times 2 ^ {1} = (1 \times 2 ^ {0} + 1 \times 2 ^ {-1} + 1 \times 2 ^ {-2} + 1 \times 2 ^ {-3}) \times 2 ^ {1}

然后得到浮点数的计算公式。

x = \pm 2 ^ {e_{x}} (1 + m_{x}) \tag{2}

$3.75_{10}$ ，就可以写为。

3.75_{10} = 2 ^ {1} (1 + 0.875)

$log_2$ 方程。

log_2(x) = e_x + log_2(1 + m_x)

$m_x \in [0,1)$ $\sigma \in [0,1)$ $\sigma = 0.0430357$ 是该算法的最优参数。

log_2(1 + m_x) \approx m_x + \sigma

该公式可以近似替换为。

log_2(x) \approx e_x + m_x + \sigma \tag{3}

然后第二步我们将float类型二进制转换为int类型数值，首先先看一下float的二进制格式。

根据浮点数编码的定义。

$S = B_{31}，0为正，1则为负$
$E = e_x + B，E \in [0,255]，B = 127，e_x \in [-127,127]$
$M = m_x \times L，L = 2 ^ {23}$

$3.75_{10} = 1.111 \times 2 ^ {1}$ 相对应的值为，可以用这个在线的浮点数转换网站验证我们的结果。

$S = 0$
$E = 1 + 127$
$M = 0.111 \times 2 ^ {23} = (0 \times 2 ^ {0} + 1 \times 2 ^ {-1} + 1 \times 2 ^ {-2} + 1 \times 2 ^ {-3}) \times 2 ^ {23} = 0.875 \times 2 ^ {23}$

于是我们可以得到如何解析浮点数二进制转换成整数的公式。

\begin{align*} I_x &= {S} \times 2 ^{31} + E \times L + M \\ &= 0 + (e_x + B) \times L + m_x \times L \\ &= L(e_x + m_x + B) \end{align*}

根据公式(3)得。

\begin{align*} I_x &= L(e_x + m_x + \sigma + B - \sigma) \\ &\approx L(log_2(x) + B - \sigma) \\ log_2({x}) &\approx \frac {I_x}{L} - (B - \sigma) \end{align*}

$I_x$ $L(log_2(x) + B - \sigma)$ 的plot曲线是很近似的。

Log_by_aliasing_to_int

带入公式(1)得。

log_2(y) = - \frac{1}{2} log_2(x) \\ \frac {I_y}{L} - (B - \sigma) \approx - \frac{1}{2} \Big(\frac {I_y}{L} - (B - \sigma) \Big) \\ I_y \approx \frac{3}{2} L(B - \sigma) - \frac{1}{2} I_x

$L, B, \sigma$ 是常量，于是我们可以得到代码里的"magic number"。

\frac{3}{2} L(B - \sigma) = 0x5f3759df \\ \frac{1}{2} I_x = i << 1


xxxxxxxxxx
i = 0x5f3759df - ( i >> 1 ); // what the fuck?

注释3: 1st iteration


xxxxxxxxxx
y = y * ( threehalfs - ( x2 * y * y ) );

Newton's method

\begin{align*} f'(x_n) &= \frac{f(x_n) - 0}{x_n - x_{n+1}} \tag{4} \\ x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \end{align*}

$y = \sqrt{x}$ 。

\begin{align*} y &= \sqrt{x} \\ y ^ 2 - x &= 0 \\ f(y) &= y ^ 2 - x \\ f'(y) &= 2y \end{align*}

带入公式(4)。

\begin{align*} y_{n+1} &= y_n - \frac{y_n ^ 2 - x}{2y_n} \\ &= \frac{y_n ^ 2 + x}{2y_n} \\ &= \frac{y_n ^ 2 + x}{2y_n} \\ &= \frac{1}{2} (y_n + \frac{x}{y_n}) \end{align*}


xxxxxxxxxx
float sqrt_impl_by_newton_method(float x)
{
    // initial guess
    float y = 0.5f * x;
    for (int i = 0; i < 5; i++)
    {
        y = 0.5f * (y + x / y);
    }
    return y;
}

$y = \frac{1}{\sqrt{x}}$ 。

\begin{align*} y &= \frac{1}{\sqrt{x}} \\ y ^ {-2} - x &= 0 \\ f(y) &= y ^ {-2} - {x} \\ f'(y) &= \frac{-2}{y^3} \end{align*}

带入公式(4)。

\begin{align*} y_{n+1} &= y_n - \frac{y_n ^ {-2} - {x}}{\frac{-2}{y_n^3}} \\ &= \frac{\frac{-2}{y_n^3} \times y_n}{\frac{-2}{y_n^3}} - \frac{y_n ^ {-2} - {x}}{\frac{-2}{y_n^3}} \\ &= \frac {-3y_{n} ^ {-2} + x} {\frac{-2}{y_n^3}} \\ &= \frac {-3y_{n} - x \times y_n ^ {3}} {-2} \\ &= y_n \times (1.5 - 0.5 \times x \times y_n ^ {2}) \end{align*}

源代码的实现。


xxxxxxxxxx
y = y * ( threehalfs - ( x2 * y * y ) );

最后

wiki有一篇很棒的文章关于这个算法实现，以及讲述科学计数法和浮点数编解码是如何相关联的一篇blog，祝阅读愉快。