Skip to content
Lex Blog

快速平方根倒数算法

本文介绍一种名为快速平方根倒数（Fast Inverse Square Root）的算法，用于快速计算浮点数的平方根倒数。

引言

快速平方根倒数（Fast Inverse Square Root）是一种相当著名的算法，其具体作者已无从考证。它因在游戏《雷神之锤 III》（Quake III）中被大量使用而闻名。

为什么需要平方根倒数而不是平方根本身呢？在计算机图形学中，许多计算涉及的是平方根倒数而非平方根。例如，向量的 normalize() 操作 1x2+y2+z2(x,y,z)\frac{1}{\sqrt{x^2 + y^2 + z^2}}(x,y,z) 就包含平方根倒数。在计算机图形学中，每一帧的渲染至少需要数万次 normalize 调用。若能加速平方根倒数的计算，渲染速度将显著提升。

该算法诞生于 20 多年前。时至今日，它只有研究价值而无实际用途，因为现代 CPU 已有更快、更精确的指令。我们直接使用 1/sqrt(y) 即可。

然而，这并不妨碍我们体会该算法的精妙之处。简要来说，它涉及：

  1. IEEE 754 标准
  2. C 语言中的未定义行为
  3. 牛顿迭代法

这里重点讨论第二步——如何推导出魔数 0x5f3759df。更详细的信息可参阅所引资料。

算法

float Q_rsqrt(float number)
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = number * 0.5F;
    y = number;

    /* 核心代码 */
    // 1. 将 float 重新解释为 int
    i = *(long *)&y; // 邪恶的浮点位级技巧
    // 2. 估计平方根倒数
    i = 0x5f3759df - (i >> 1); // 这到底是什么鬼？(....为什么不是 0x69696969?)
    // 将 int 重新解释为 float
    y = *(float *)&i;
    // 3. 牛顿法
    y = y * (threehalfs - (x2 * y * y)); // 第 1 次迭代
    // y = y * (threehalfs - (x2 * y * y)); // 第 2 次迭代，可以删除

    return y;
}

该算法大致由三部分组成：

  • 用 reinterpret_cast 在 float 与 int 之间转换，因为任何编译器都会禁止对 float 类型进行位运算。
  • 利用魔数 0x5f3759df 估计平方根倒数（后文说明）。
  • 牛顿法的一次迭代（后文说明）。

估计平方根倒数

在 IEEE 754 标准中，一个 32 位浮点数包含 1 位符号位、8 位指数位 EE 与 23 位尾数位 MM。任何浮点数都可以表示为：(1.M)×2E−127=(1+M223)×2E−127,M∈[0,223),E∈[0,256)(1.M) \times 2^{E - 127} = (1 + \frac{M}{2^{23}}) \times 2^{E - 127}, M \in[0, 2^{23} ), E \in [0,256)。

为了计算 yy 的平方根倒数 1y=y−12\frac{1}{\sqrt{y}} = y^{-\frac{1}{2}}，yy 的 IEEE 754 形式为 y=(1+My223)×2Ey−127y = (1 + \frac{M_y}{2^{23}}) \times 2^{E_y - 127}。

为简单起见，对 yy 取以 2 为底的对数：

log⁡y=log⁡((1+My)×2Ey−127)=log⁡(1+My223)+Ey−127\begin{align} \log y &= \log{((1 + M_y) \times 2^{E_y - 127})} \\ &= \log{(1+\frac{M_y}{2^{23}})} + E_y - 127 \end{align}

显然 My223\frac{M_y}{2^{23}} 位于区间 [0,1][0,1] 内。在这个区间上，log⁡(1+x)≈x\log(1+x) \approx x。若加上一个常数 μ\mu，近似效果会更好：

利用这个近似，可以进一步化简：

log⁡y=log⁡(1+My223)+Ey−127≈My223+μ+Ey−127=1223(My+Ey×223)+μ−127\begin{align} \log y &= \log{(1+\frac{M_y}{2^{23}})} + E_y - 127 \\ &\approx \frac{M_y}{2^{23}} + \mu + E_y - 127 \\ &= \frac{1}{2^{23}}(M_y + E_y \times 2^{23}) + \mu - 127 \end{align}

现在，对 y−12y^{-\frac{1}{2}} 取对数：

log⁡y−12=−12log⁡y≈−12(1223(My+Ey×223)+μ−127)\log y^{-\frac{1}{2}} = -\frac{1}{2} \log y \approx -\frac{1}{2}\left( \frac{1}{2^{23}}(M_y + E_y \times 2^{23}) + \mu - 127 \right)

假设 y−12y^{-\frac{1}{2}} 的精确值为 AA，其指数与尾数分别为 EAE_A 与 MAM_A，可以对 AA 取对数：

log⁡A≈1223(MA+EA×223)+μ−127\log A \approx \frac{1}{2^{23}}(M_A + E_A \times 2^{23}) + \mu - 127

令两个对数表达式相等，得到：

MA+EA×223=223×32(127−μ)−12(My+Ey×223)\begin{align} M_A + E_A \times 2^{23} &= 2^{23} \times \frac{3}{2}(127 - \mu) - \frac{1}{2}(M_y + E_y \times 2^{23}) \end{align}

当 μ\mu 取 0.04504618757916870117560.0450461875791687011756 时，第一项变为十六进制的 1597463011（0x5F3759e3）。若把第二项的 12\frac{1}{2} 换成右移操作，结果与代码中的那行惊人地接近：

i = 0x5f3759df - (i >> 1);

魔数略有不同，但只要 μ\mu 取值合适，就能得到完全相同的数。

牛顿迭代法

牛顿迭代法的更新公式为：

xt+1=xt−f(xt)f′(xt)x_{t+1} = x_t - \frac{f(x_t)}{f^\prime(x_t)}

牛顿迭代法的前提是初始点应靠近函数的零点。由于 i = 0x5f3759df - (i >> 1); 得到的估计值误差很小，我们可以认为该前提已满足。

现在的问题是：给定 yy，需要求 xx 使得 1y=x\frac{1}{\sqrt{y}} = x，即 1x2=y\frac{1}{x^2} = y。整理得 f(x)=1x2−y=0f(x) = \frac{1}{x^2} - y = 0。对 f(x)f(x) 应用牛顿迭代，得到迭代公式：

xt+1=xt−y2xt3x_{t+1} = x_t - \frac{y}{2}x_t^3

翻译成代码即为：

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

参考资料

1 Is fast inverse square root still used?

2 Fast Inverse Square Root — A Quake III Algorithm

3 Relevant Thesis