快速平方根倒数：0x5f3759df 这个魔法数是怎么来的

在《Quake III Arena》的源代码里，有一段非常出名的函数。它的目标很简单：对于一个正浮点数 \(x\)，快速近似计算

代码却一点也不普通：

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

    return y;
}

最令人困惑的显然是这一行：

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

为什么把一个浮点数的二进制位当成整数，右移一位，再从一个十六进制常数中减掉，就能得到平方根倒数的近似值？为什么后面再乘一次

y * (1.5F - 0.5F * number * y * y)

结果就会突然变得很准？

整个算法可以拆成两步：

1. 利用 IEEE 754 浮点数的结构，以极低成本构造一个不错的初始值；
2. 用一次牛顿迭代快速修正初始值。

所谓“魔法”，其实是数据表示、对数近似和数值迭代恰好在这里接到了一起。

为什么要计算平方根的倒数

在三维图形中，经常需要把一个向量归一化。对向量

单位向量为

这里既有平方根，又有除法。在早期硬件上，这两类运算都比较昂贵，而向量归一化又会在光照、反射和物理计算中频繁出现。因此，如果能快速得到

再分别与三个分量相乘，就能节省大量计算。

这也是快速平方根倒数算法曾经极具吸引力的背景。不过，理解它并不需要了解游戏引擎，只需要先看看一个 float 在内存中是什么样子。

IEEE 754 单精度浮点数

常见的 IEEE 754 单精度浮点数，也就是 binary32，使用 32 个二进制位：

可以把它写成：

31          30              23 22                         0
+--------------+---------------+---------------------------+
| sign (1 bit) | exponent (8)  | fraction / mantissa (23) |
+--------------+---------------+---------------------------+

本文只推导正的正规数，因此符号位 \(s=0\)，指数满足 \(0<E<255\)。此时浮点数的实际数值是

其中 127 是指数偏置。尾数字段虽然只有 23 位，但正规数的最高位有一个没有显式存储的 \(1\)，所以有效数字写成

为了让后面的式子更清楚，记

于是：

同一串二进制位，也可以看成整数

现在暂时忘掉这 32 位表示的是 float，直接把同一串位解释成一个无符号整数。由于符号位为 0，这个整数可以写成

两边除以 \(2^{23}\)：

这里一定要区分两个量：

• \(x\) 是浮点数所表示的实际数值；
• \(I_x\) 是同一串 32 位二进制按整数解释后的数值。

原始代码中的指针转换，做的就是在这两种解释之间切换，而不是把浮点数进行普通的数值类型转换。

为什么位模式近似于对数

对

取以 2 为底的对数：

在 \(0\le m_x<1\) 的区间里，曲线 \(\log_2(1+m_x)\) 与直线 \(m_x\) 并不完全相等，但两者比较接近：它们在 \(m_x=0\) 时相等，并在 \(m_x\to1\) 时趋向同一个值 1。因此可以作近似：

于是：

这就是位运算能够发挥作用的核心：

对于正的 IEEE 754 正规浮点数，其整数位模式 \(I_x\) 可以近似看成 \(\log_2 x\) 的线性编码。

换一种直观说法，如果我们定义一种“理想对数格式”

那么它与真实浮点值

的区别只在于用 \(2^{m_x}\) 代替了 \(1+m_x\)。二者不相等，但足够接近，可以先用前者快速求一个粗略答案，再进行修正。

从 \(x^{-1/2}\) 推出“常数减去右移一位”

我们要求的是

取对数后，乘方变成了乘法：

分别用位模式的线性近似表示 \(x\) 和 \(y\)：

整理得到：

也就是

而

所以最朴素的位级初值公式是

I_y = 0x5f400000 - (I_x >> 1);

现在那行神秘代码已经有了轮廓：

• I_x >> 1 相当于对整数位模式除以 2，对应对数关系中的系数 \(-\frac12\)；
• 0x5f400000 负责补回指数偏置 127，并把结果移到正确的数量级；
• 得到的整数位模式重新解释为 float 后，就是 \(1/\sqrt{x}\) 的初始近似。

如果浮点数真的是一种完全对数化的编码，0x5f400000 就足以给出精确答案。但 IEEE 754 浮点数在相邻的 2 的幂之间是线性的，而真正的对数格式是指数曲线，所以仍然存在系统误差。

为什么代码用的是 0x5f3759df

0x5f400000 来自最简洁的理论模型，经典代码却使用：

0x5f3759df

原因是对数近似

并不精确。等价地说，在 \(0<m<1\) 时，真实浮点有效数字 \(1+m\) 与理想对数格式中的 \(2^m\) 有偏差。直接使用 0x5f400000 得到的初值会整体偏高，因此将常数稍微调小，可以让正负误差更加平衡，从而降低最大误差。

所以 0x5f3759df 并不是某个严格恒等式自然算出的唯一答案。更准确的理解是：

它是“常数减去一半位模式”这一族近似中的一个误差表现很好的参数。

而且“最好”还取决于优化目标：

• 如果只考察位运算得到的初始值，某个常数可能最好；
• 如果考察一次牛顿迭代后的误差，最佳常数可能不同；
• 如果考察两次迭代或使用另一种误差指标，答案还可能变化。

Chris Lomont 对这一问题做过系统分析。他得到：只优化初始近似时，0x5f37642f 的最大误差略小；把牛顿迭代也计算在内时，数值搜索得到的 0x5f375a86 又略优于经典常数。不过，这些数都非常接近，经典的 0x5f3759df 本身已经相当优秀。

这也解释了为什么仅靠一个简单推导，通常只能得到接近的 0x5f400000 或其他相邻常数，而不会自动、唯一地落在 0x5f3759df 上：理论给出了常数的结构和大致范围，具体低位则来自对误差的进一步优化。

牛顿迭代：把粗略答案迅速修准

位运算只负责给出一个初始近似。接下来要用牛顿迭代把它修正。

牛顿法的一般形式是

这里输入 \(x\) 是固定的，我们想找到

它等价于

因此可以构造一个以 \(y\) 为未知数的函数：

目标答案正是 \(g(y)=0\) 的正根。求导得到

代入牛顿公式：

这正是代码中的：

y = y * (1.5F - x2 * y * y);

其中 x2 = number * 0.5F，所以 x2 就是 \(x/2\)。

牛顿法在初值已经比较接近答案时收敛很快。Lomont 的测试中，经典常数产生的初始值最大相对误差约为 \(3.44\%\)；一次牛顿迭代后降到约 \(0.175\%\)；第二次迭代后又下降了几个数量级。这就是算法先花极低成本“猜一个差不多”，再用一次迭代迅速逼近答案的原因。

用现代 C++ 安全地实现

经典代码写于特定的编译器和硬件环境，有两个可移植性问题：

1. C/C++ 并不保证 long 一定是 32 位；
2. 把 float* 强制转换成 long* 后解引用，会触犯现代 C++ 的严格别名规则，行为未定义。

C++20 提供了 std::bit_cast，可以在保持位模式不变的前提下安全地切换解释方式。配合固定宽度整数，代码可以写成：

#include <bit>
#include <cstdint>

float fast_inverse_sqrt(float x) {
    const float half_x = 0.5F * x;
    std::uint32_t bits = std::bit_cast<std::uint32_t>(x);
    bits = 0x5f3759dfU - (bits >> 1);

    float y = std::bit_cast<float>(bits);
    y = y * (1.5F - half_x * y * y);
    return y;
}

这个版本避免了指针类型双关，但算法本身仍然依赖两个前提：float 使用 32 位 IEEE 754 binary32 表示，并且输入处在算法所假设的范围内。

输入范围与特殊值

前面的推导一直假设 \(x\) 是正的正规浮点数。以下输入不能直接套用这套推导：

• \(x=0\)：数学上的 \(1/\sqrt{x}\) 发散；
• \(x<0\)：在实数范围内没有平方根；
• 正负无穷与 NaN：指数位具有特殊含义；
• 次正规数：没有正规数所使用的隐含最高位 \(1\)，指数公式也不同。

因此，上面的函数更适合作为算法原理的现代化表达，而不是一个对所有 float 输入都定义良好的标准库替代品。如果用于真实项目，应当明确前置条件，或者在外层处理特殊值。

今天还应该使用它吗

不能脱离平台直接断言它一定比

1.0F / std::sqrt(x)

快多少。今天的处理器可能提供硬件平方根或倒数平方根指令，编译器也可能在允许快速数学优化时自动生成合适的 SIMD 指令。另一方面，标准库版本、编译选项、目标精度以及数据是否便于向量化，都会影响最终结果。

如果要在现代工程中使用类似技巧，至少要同时测量两件事：

1. 在目标硬件和真实数据上的性能；
2. 最大误差、平均误差以及特殊输入的行为是否满足需求。

从实用角度看，它已经不再是一个应该无条件复制的性能秘方；从思想上看，它依然是一段极其漂亮的代码。

总结

快速平方根倒数可以概括成三步：

1. IEEE 754 浮点数的整数位模式近似线性编码了它的二进制对数；
2. 在对数坐标中，\(x^{-1/2}\) 只需要乘以 \(-1/2\)，于是变成“魔法常数减去右移一位”；
3. 位运算给出初值后，再用一次牛顿迭代快速降低误差。

因此，0x5f3759df 真正精彩的地方不是它看起来神秘，而是它把三个不同层次的知识压缩进了几行代码：底层的数据表示、中间的数学近似，以及上层的数值迭代。

需要注意的是，这是一种根据代码结构进行的数学解释。它说明了为什么这一族算法能够成立，以及常数应当落在什么范围；但仅凭现有资料，不能反推出原作者当年究竟是完全通过推导得到常数，还是结合了试验和搜索。

参考资料

• id Software：Quake III Arena 中的 Q_rsqrt 源码
• Oracle Numerical Computation Guide：IEEE 单精度格式
• Chris Lomont：Fast Inverse Square Root
• Tim Hutt：Fast Inverse Square Root
• Better Explained：Understanding Quake’s Fast Inverse Square Root
• cppreference：std::bit_cast