- 0x5f3759df 是快速平方根倒数算法的核心,它通过位运算直接给出 1/√x 的初始近似值。
- 配合牛顿迭代法,只需 1~2 次迭代就能达到极高精度,整体速度超传统 sqrt 。
- 这种“位级黑科技”是当年程序员在硬件受限下的极致优化,至今仍是计算机科学中的经典案例。
一、神奇的 0x5f3759df
在快速平方根倒数算法中, 0x5f3759df 是那个让无数程序员拍案叫绝的魔法数字(Magic Number)。它能在没有浮点运算单元(FPU)的年代,用纯整数位运算就给出平方根倒数的一个极佳初始近似值,再配合牛顿迭代法快速收敛到精确结果,整体速度比传统 sqrt 快数倍。
二、算法核心代码解析
float Q_rsqrt( float number ) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // 邪恶的浮点位级黑科技
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 第一次迭代
// y = y * ( threehalfs - ( x2 * y * y ) ); // 第二次迭代,可选
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}
这段代码的关键思路:
1. 位级黑科技:把浮点数的二进制位直接当成整数来操作,这是实现近似的基础。
2. 魔法数字计算: 0x5f3759df - (i >> 1) 一步就给出了 1/√number 的初始估计值。
3. 牛顿迭代法:用 y = y * (1.5 - (x2 * y * y)) 迭代修正,一次迭代就能把误差降到千分之二以内。
三、浮点数格式与初始估计原理
要理解魔法数字的来源,先看 32 位单精度浮点数的结构:
- 符号位 s(1 位):0 为正,1 为负
- 指数位 E(8 位):存储 E - 127 (偏移量 127)
- 尾数位 M(23 位):存储小数部分,实际值为 1 + M/2²³
浮点数的数值公式:
x = (-1)^s \left(1 + \frac{M}{2^{23}}\right) 2^{E-127}
我们想求 y = \frac{1}{\sqrt{x}},两边取二进制对数并近似,最终可以推导出:
- 把浮点数 x 的整数表示 i 右移一位,再用魔法数字 0x5f3759df 去减,就能得到 y 的整数近似表示。
- 这个魔法数字是通过数学推导和大量实验找到的,能让初始估计的相对误差最小。
四、魔法数字的求解:三分搜索
魔法数字不是凭空来的,而是通过三分搜索在合理区间内找到的最优值。
if __name__ == "__main__":
mant = np.arange(0, 1 << 23, dtype=np.uint32)
# 给一个合理的搜索区间,围绕 0x5f3759df 一带
lo = 0x5f000000
hi = 0x5f900000
best_magic, best_err = ternary_search_magic(lo, hi, mant)
print(f"best_magic = 0x{best_magic:08x}")
print("最小最大相对误差 =", best_err)
实验结果对比(以 32 位为例):
来源 值 初始估计相对误差 一次迭代相对误差
数学证明值 0x5f37642f 3.421282e-02 1.775889e-03
三分查找 0x5f375a87 3.436540e-02 1.751287e-03
Magic Number 0x5f3759df 3.437577e-02 1.752338e-03
64位数学证明值 0x5fe6ec85e7de30da 3.421281e-02 -
可以看到,0x5f3759df 与理论最优值几乎一致,是工程上的最优选择。
五、实验效果:64位对比 Magic Number
在 64 位双精度浮点数中,也有对应的魔法数字 0x5fe6ec85e7de30da 。它同样遵循位运算近似 + 牛顿迭代的思路,只是浮点数的指数位(11 位)和尾数位(52 位)长度不同,魔法数字的取值也随之调整