黑魔法不是魔法：开方为何变成减法

1999 年 Quake III Arena 的源码里，藏着一行让全世界程序员愣了二十多年的注释。它出现在一个叫 Q_rsqrt 的函数中，紧跟在一行奇怪的整数运算后面：

i  = * ( long * ) &y;                // evil floating point bit level hacking
i  = 0x5f3759df - ( i >> 1 );        // what the fuck?

// what the fuck?——连写下它的人都承认自己没完全搞懂这是怎么回事。这段代码要做的事极其普通：算 1/√x，即一个数的倒数平方根（inverse square root）。可它既没有除法，也没有 sqrt，只用了一次整数右移、一次减法、一次乘法，就蹦出了一个精度足够做 3D 光照的答案，比当年标准的浮点除法快好几倍。

ThisPrimeagen（ThePrimeTime 频道）的反应视频《Prime React: Fast Inverse Square Root — A Quake III Algorithm》之所以能拿下 25 万+播放，正是因为它戳中了每个程序员的同一个疑问：一串来历不明的十六进制常数，凭什么替你把开方算出来了？ 这个问题的答案，是理解"数据表示"这堂课最精彩的一份教材。下面我把这层"黑魔法"一层层剥开给你看——你会发现，它根本不是魔法，而是浮点数偷偷藏好的一本对数表。

视频背景卡片

一段源码引出的悬案：它到底在干什么

项目	内容
视频标题	Prime React: Fast Inverse Square Root — A Quake III Algorithm
频道	ThePrimeTime（ThePrimeagen，百万订阅级编程博主）
时长	约 25 分钟（reaction + 讲解）
播放量	约 25.6 万
发布日期	2023-02-03
核心主题	拆解 Quake III 的 Fast Inverse Square Root 算法与 `0x5f3759df` 常数

先把整个函数摆出来（来自 id Software 在 GitHub 上的开源仓库 code/game/q_math.c 第 552 行）：

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;
}

读法是关键：&y 取浮点变量 y 的地址，(long*)&y 把它"骗"成整数指针，再解引用——于是一个 32 位浮点数被原封不动地当成 32 位整数 i 来读，一位都不动。这就是注释里说的 "bit level hacking"：不经过任何类型转换函数，直接换一副眼镜看同一串比特。接着 i >> 1 把这个整数右移一位（除以 2），用常数 0x5f3759df 减掉它，再把结果按浮点重新解释回去。最后一行乘法是 Newton 迭代的"擦亮"步骤。

问题来了：把浮点当整数右移再相减，凭什么等价于开方？这串魔法常数又是哪来的？要回答它，必须先弄懂 IEEE 754 浮点数在内存里到底长什么样。这也是为什么这段代码被无数人奉为"理解计算机数值表示"的最佳入口——它逼你看穿表象。

第一根支柱：浮点数的"位"，偷偷就是它的对数

一个正的 IEEE 754 单精度浮点数 x，在内存里是 1 位符号、8 位指数 E、23 位尾数 M 拼成的 32 位，其值为 x = (1 + M/2²³) · 2^(E−127)。把它当作一个整数 I 来读，则 I = E·2²³ + M。现在做个神来之笔：对 x 取以 2 为底的对数。

log₂ x = (E − 127) + log₂(1 + M/2²³)

而 log₂(1+t) 在 t ∈ [0,1) 上近似一条直线，约等于 t + σ（σ ≈ 0.043，是让误差最小的一个小常数）。代入并整理后，你会得到一条极其漂亮的关系：

I  ≈  L · log₂ x + B       （其中 L = 2²³，B = 2²³·(127 − σ)）

也就是说，把浮点数当成整数读，这个整数值几乎就是 log₂ x 的线性函数。浮点数在自己的位模式里，悄悄存着一份自己的对数。这是一切魔法的基础：既然 I 线性正比于 log₂ x，那么对 I 做的算术运算，就等价于对 log₂ x 做同样的运算。位运算和对数运算，在这里是同一件事的两副面孔。

第二根支柱：减一个常数，等于把对数乘以 −½

我们想算的是 1/√x = x^(−½)。两边取对数：log₂(x^(−½)) = −½ · log₂ x。目标浮点数 y = 1/√x 的整数位模式 I_y，按上一节的线性关系应为：

I_y  ≈  L · (−½ log₂ x) + B
     =  −½ (I_x − B) + B
     =  −½ · I_x  +  (3/2)·B

看清楚了吗？右边的 −½ · I_x，正是源码里的 i >> 1（右移一位 = 除以 2）；而 (3/2)·B，正是那串魔法常数 0x5f3759df。所以这行看起来像胡来的整数减法：

i = 0x5f3759df - ( i >> 1 )

在数学上等价于"把 log₂ x 乘以 −½"——也就是在算 1/√x。 常数不是魔法，它就是线性关系里的那个截距项 (3/2)·B；右移不是魔法，它就是把指数（也就是对数）减半。整段"黑魔法"翻译成人话只有一句：浮点数已经替你存好了对数，开方不过是把对数除以二，再取个负号。 难怪原作者写下后忍不住骂了一句。

那串常数不是乱猜的：Lomont 把它推了回来

民间一直传说 0x5f3759df 是某位天才"试出来的"。2003 年，普渡大学的 Chris Lomont 写了一篇严谨的论文《Fast Inverse Square Root》，把这事彻底说清：常数完全可由上面的线性关系解析推出，而非瞎猜。关键在那个偏置常数 σ——选 σ 让"最大相对误差最小"（minimax 准则），就能反推出最优常数。

Lomont 算出经典常数 0x5f3759df 给出的初值最大相对误差约为 0.175%（1.75×10⁻³），并且他还找到了一个略优的替代 0x5f375a86，把误差再压低一丁点。这等于给"魔法"盖了章：它不仅可解释，还能用中学水平的近似加上一点点优化理论从头推导。所谓黑魔法，其实是被人讲清楚之后就索然无味的初等数学。这段历史也是想提醒你：看到一个"魔法常数"，第一反应不应该是敬畏，而是追问它对应的数学结构是什么。

Newton：真正把精度擦亮的，是最后那行乘法

常数那一步给出的只是"种子"，离精确还差得远。真正把答案打磨到位的，是最后一行的 Newton–Raphson 迭代：

y = y * ( threehalfs - ( x2 * y * y ) );   // x2 = x/2

令 g(y) = 1/y² − x，它的根就是 y = 1/√x。对它套牛顿法 y_new = y − g(y)/g'(y)，化简后正好是 y·(1.5 − (x/2)·y²)，与代码完全一致。牛顿法是二次收敛的：误差每迭代一次大致"平方一次"。所以初值那 0.175% 的相对误差，经过一次迭代就塌缩到约 10⁻⁶ 量级；源码里注释掉的第二行，是给你想再迭代一次、把精度推到约 10⁻¹¹ 用的。

这就是这段代码的精妙分工：位运算负责"几乎免费的、误差千分之一的初值"，牛顿法负责"把误差平方掉"。 缺了位运算那步，牛顿法从一个烂起点要迭代很多次；缺了牛顿法，位运算的初值又不够精确。两者缺一不可，组合起来才又快又准。

谁写的？连 Carmack 都说"不是我"

这串代码的署名之争本身就像侦探小说。因为 Quake III 的灵魂人物是 John Carmack，多年来所有人都默认是他写的。但 Carmack 本人公开否认，说它早就在代码库里了。后来 Beyond3D 论坛上一场持续数年的考据（"Origin of Quake3's Fast InvSqrt()"，Rys Sommefeldt 主导）一路追到了 Gary Tarolli（Nvidia），再追到 Greg Walsh（1980 年代末 Ardent Computer 的 Titan 图形工作站），其思想根源甚至能上溯到数值计算大师 William Kahan。换句话说，这招比 Quake III 老了十来年，是在大型图形工作站上为了给光照归一化向量而发明的。

这段历史的价值在于：软件世界里那些看似"个人天才闪光"的东西，往往是一条绵延几十年的工程传统。把它当成孤立的传奇，你会觉得可遇不可求；把它放回传统里，你会发现这是一类问题（"如何在受限硬件上用近似换速度"）反复被打磨出的标准答案。这也解释了为什么它值得每个程序员都看一遍——它代表的不只是一个技巧，而是一种思维方式。

今天还需要它吗：硬件早就替你做了

诚实地说：在现代桌面/移动 CPU 上，这段 C 代码基本是"教学遗产"了。Intel SSE 的 rsqrtss（约 3–4 周期、12 位精度）、AVX-512 的 vrsqrt14ss（14 位精度）、ARM 的 vrsqrtes_f32（约 2–3 周期）都在硬件层面提供等价甚至更优的指令。当年那段位运算省下的几十个浮点周期，今天硬件用一个周期就还给你了。所以别在 2026 年的项目里照抄它"优化"——大概率比编译器+硬件指令更慢、更不精确。

但它的思想血脉从未断过。"浮点的位 ≈ 它的对数" 这个洞察，正是当代一大类技术的基础：对数数制（Logarithmic Number System）、低功耗神经网络训练的 LogNet、OCP 微缩放格式（Microscaling/MX，把指数分组共享以压缩 LLM 权重）……甚至当下 4-bit 量化大模型的底层直觉——"动态度量藏在指数里"——和 Quake 这行代码是同一个祖先。学会它，你抓到的不是一段过时代码，而是一把理解现代 AI 数值压缩的钥匙。

知识结构图

flowchart TD
    A["目标: y = 1/√x = x^(−½)"] --> B["取对数: log₂y = −½·log₂x"]
    B --> C["洞察: 浮点位模式 I ≈ L·log₂x + B<br/>(IEEE 754 偷藏对数)"]

    C --> D1["i >> 1<br/>右移一位 = log₂x 除以 2"]
    C --> D2["0x5f3759df = (3/2)·B<br/>截距项, 可解析推出"]

    D1 --> E["i = 0x5f3759df − (i>>1)<br/>== 把 log₂x 乘以 −½"]
    D2 --> E
    E --> F["重解释为浮点 → 初值<br/>最大相对误差 ≈ 0.175%"]

    F --> G["Newton 迭代一次<br/>y·(1.5 − x/2·y²)<br/>误差平方收敛 → ≈10⁻⁶"]
    G --> H["可用答案 (3D 光照归一化)"]

    style C fill:#fff4e6,stroke:#e8a33d
    style E fill:#ffe8d6,stroke:#e8590c
    style G fill:#e6f4ff,stroke:#339af0

实战启示

第一，别被"魔法常数"吓倒。 下次你在代码或论文里看到一个奇怪常数，先假设它能被解释。Quake 这串 0x5f3759df 看起来像玄学，本质只是线性近似的截距 (3/2)·B。把"这是不是某个近似公式的系数"列为第一假设，你会少很多迷信。

第二，吃透你的数据表示。 这段代码的全部威力来自"知道 IEEE 754 的位模式近似对数"。同理：知道 UTF-8 是自同步的、知道 Git 存的是快照而非差异、知道 LSM 树从不真的修改——这些"表示层的真相"往往比算法本身更值钱。真正的高性能优化，多半发生在"换个视角看同一份数据"的那一刻，而不是在循环里抠一两条指令。

第三，先粗后精是个通用范式。 位运算出初值 + 牛顿法擦亮，是"快速近似 + 局部精修"的经典组合。这个范式今天依然到处都是：哈希表用廉价哈希先定位再精确比较、网络用 CRC 先校验再重传、甚至 ML 推理用低精度算粗值再高精度校正。学会"用廉价手段拿到够好的初值，再用昂贵手段精修"，你就掌握了性能工程的一半。

第四，别照抄它"优化"。 2026 年的桌面/移动平台有硬件 rsqrt 指令，软件版只会更慢更差。它在嵌入式、没有 FPU 的 MCU、或教学场景里仍有价值；在主流平台，请直接用 rsqrt 内建函数或 std::rsqrt（若有）。判断一段"经典优化"是否还有效，要看它优化的瓶颈今天是否还存在——浮点除法周期数的瓶颈，二十年前就消失了。

延伸学习路径

一句话总结

0x5f3759df 不是魔法常数，它是 IEEE 754 位模式里那本隐藏对数表的截距——开方在浮点的世界里，真的就只是一次减法和一次右移。