卡马克快速平方根倒数算法

背景

对于许多人而言，求解平方根倒数是一个十分简单的事情，我们只需要调用C语言中头文件便可以实现

#include<stdio.h>
#include<math.h>

int main(void){
	float x = 1 / sqrt(2);
	printf("%lf",x);
	return 0;
}

但是，有没有一种方法可以更加快速地运算出结果呢

雷神之锤

在2005年，id Software 发布了 Quake-III Arena（雷神之锤3）游戏的源代码，在那份源代码中，游戏粉丝们发现了一个算法十分巧妙，这个算法的唯一用途就是计算平方根的倒数
$$
f(x)=\frac1{\sqrt{x}}
$$

这款游戏在上世纪九十年代风靡全球，即使当时的计算机配置低，但也能流畅地运行这款画面和内容都不错的游戏。这其中最重要的原因就是它的 3D 引擎的开发者 John Carmack（约翰-卡马克）

在游戏里，要想实现物理效果、光影效果，反射效果的话，要对使用的向量进行单位化，因而计算平方根倒数是一项重要且关键的步骤。
$$
z\frac1{\sqrt{x^2+y^2+z^2}}
$$
可以看见，在实现单位化的过程中，倒数平方根是一个重要的计算量。

事实上早在90年代初DOS时代，只要能在PC上搞个小动画都能让人惊叹一番的时候，John Carmack就推出了石破天惊的Castle Wolfstein, 然后再接再励，doom, doom-II, Quake…每次都把3D技术推到极致。他的3D引擎代码极度高效，几乎是在压榨PC机的每条运算指令。当初微软的Direct3D也得听取他的意见，修改了不少API。

而在源代码中，越是底层的函数调用的就越频繁，就越是一个游戏引擎优秀的关键，在 Quake-III Arena 的源代码文件 game/code/q_math.c 里就有这样一段神奇的平方根倒数的算法

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
    #ifndef Q3_VM
    #ifdef __linux__
         assert( !isnan(y) ); // bk010122 - FPE?
    #endif
    #endif
    return y;
}

经过测试比较，一般而言卡马克快速平方根算法速度是普通算法 sqrt 的三到四倍

原理分析

首先你会注意到一个奇怪的数字 0x5f3759df ，这个数字是怎么样得出来的呢，这一步又有什么用途呢？

而且，作者在注释中还写了很有趣的一句话 what the fuck?

二进制数的分析

观察代码的头三行可以发现：

1、作者首先定义了一个 long 类型的数 i ，也就是说，我们现在有32个二进制位来表示数字

2、在下一行作者声明了两个小数 float ，也就是说，我们现在有32个二进制位来表示小数。

小数的二进制表示

一般来说，类似于整数的二进制表示方法，我们会采用将32位二进制数分割的办法来表示小数，例如：00000000 00000001.10000000 00000000 表示 1.5。

但这样会有一个问题，原本可以表示 $2^{31}$ 大小的数，现在只能表示 $2^{15}$ 大小的数了，这造成了严重的资源浪费。所以为了更好的表示小数，计算机科学家们借鉴了科学计数法的表示方法，制定了 IEEE-754 标准

IEEE-754 标准

对于32位的单精度浮点数（float）存储，按照IEEE-754的标准总共分为三个部分：

• 符号位 S ：1 位，0 表示正数，1表示负数
• 指数位 E ：8 位，127 偏移量
• 尾数位 M ：23 位

示例

下面我们来看一个简答的例子，了解一下计算机是如何存储浮点数的

对于二进制数 0.101 相当于十进制数 0.625

我们可以类比科学计数法的思想来表示 0.625
$$
(-1)^{0}\times2^{-1}\times1.01
$$
其中 $(-1)^{0}$ 是符号位；$2^{-1}$ 是指数位，底数都是 $2$；$1.01$ 是尾数位

对于真正的二进制表示 0.625 可得 0 01111110 01000000000000000000000

指数位分析

首先我们分析指数位，他用二进制 01111110 来表示 2 的指数 -1，那这其中的原因是什么呢？

01111110 表示 126 ，而 8 位的二进制总共可以表示的大小范围为 255 (11111111)

查询进制转化表可知

-2 对应 125

-1 对应 126

1 对应 127（0 用其他方式表示）

2 对应 128

也就是这样一种对应方式
$$
2^x = 2^{127 + x}
$$

尾数位分析

如果对科学计数法熟悉的话，可以发现最后的数值大小范围始终在 1 ~ 2 之间，所以我们可以将32位全部用来表示小数部分，最后在计算机的底层处理中再加入1就可以了。

例如，1.01 中只需要表示0.01就可以了，所以就直接用二进制表示 01000000000000000000000

练习

大家可以思考一下，对于100.25这个数，该如何用IEEE754标准来表示呢？

1. 首先我们要讲他转化成二进制数来表达 1100100.01
2. 然后可以将二进制数写成 $(-1)^{0}\times2^{6}\times1.10010001$
3. 接着确定指数位：$127+6=133$ ，133 又可以用二进制写成 10000101
4. 最后来确定尾数，直接应用2中最后的小数部分，10010001 00000 00000 00000

综上，100.25 这个数用IEEE标准的二进制表示为 0 10000101 10010001 00000 00000 00000

十六进制数的分析

如果你已经按照上诉步骤学会了IEEE754标准里的单精度浮点数的二进制表示方法，相信你可以类推到双精度浮点数（double）上，也可以将二进制的表示方法推广到十六进制的表示方法上——四位为一组进行转换

求初始近似值

卡马克快速平方根算法的核心在于牛顿迭代法的使用，但是牛顿迭代法如果初始的精度值不够，就必须使用多次的牛顿迭代法来达到更大的精度。但是，如果我们能先求出来一个精度还不错的初始值，再对它使用一次牛顿迭代法，就可以一次性得到更大的精度。所以，我们这里的关键是求出一个近似的平方根倒数初始值。

公式推导

根据前面的IEEE754标准的浮点数表示方法，我们可以将浮点数 y 写成下列公式来表示
$$
y=(1+\frac M{2^{23}})\cdot2^{E-127}
$$
M ：尾数位；E ：指数位

那么，在计算机中怎么样可以快速地求出 $y^{-\frac{1}{2}}$ 呢？我们可以采用取对数的方式来进行一些神奇的变换！
$$
\begin{aligned}&\log_2\left[\left(1+\frac{M}{2^{23}}\right)\cdot2^{E-127}\right]\&=\log_2\left(1+\frac M{2^{23}}\right)+\log_22^{E-127}\&=\log_2\left(1+\frac M{2^{23}}\right)+E-127\&\approx\frac M{2^{23}}+E-127\end{aligned}
$$
通过这样的一系列变化，我们便可以很轻松的得到 $\log_2{y}$ 的近似值，从而得到 $-\frac{1}{2}\log_2{y}$ 的值，进而得到 $y^{-\frac{1}{2}}$ 的初始近似值

进一步，我们还可以化简公式使得更接近于计算机底层的浮点数表示方法
$$
\frac1{2^{23}}(2^{23}\cdot E+M)-127
$$
其中， $2^{23}\cdot E+M$ 就是电脑里存储的浮点数 y 的二进制格式

最后，我们使用 A 来代替初始近似值的二进制表示，Y 来代替原始值的二进制表示，可以得到这样一个等式
$$
\frac{1}{2^{23}}\cdot A-127=-\frac{1}{2}\left(\frac{1}{2^{23}}\cdot Y-127\right)\A=381\cdot2^{22}-\frac{1}{2}Y
$$

代码解释

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;

第一步

我们定义了一个长整数类型的变量 i ，接着让它来表示单精度浮点数 number 的二进制格式（因为在底层它们的大小都是32位）。

在这里使用的是这样一个小技巧 i = * ( long * ) &y; 如果直接使用 i = ( long ) y 这将会改变底层的二进制格式。

第二步

套用公式 $A=381\cdot2^{22}-\frac{1}{2}Y$ 来计算初步近似值，但是如果你自己仔细观察可以发现：$381\cdot2^{22}$ 这个数的十六进制表达式为 5F40 0000 ，并非公式里的值。

这是因为我们在进行 log 运算的进行了粗糙的近似，如果使用 0x5f3759df 可以弥补这样的精度损失

第三步

最后再将 i 的底层二进制格式转换为 y 的单精度浮点数的表示

牛顿法迭代

公式

$$
x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime}(x_n)}
$$

$$
x_{n+1}=x_n-\frac{x_n^{-2}-x}{-2x_n^{-3}}
$$

$$
x_{n+1}=\frac{3}{2}x_n-\frac{1}{2}xx_n^3
$$

其中，$x$ 代表输入的变量 number ，$x_n$ 代表初始的近似值 $y$ ，$x_{n+1}$ 代表一次牛顿迭代之后的近似值

通过牛顿迭代法，我们可以初步逼近最后的高精度近似值

代码实现

   const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
   // y   = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

可以看到，原来的代码中还有第二次牛顿迭代的计算，但是通过实验验证了一次牛顿迭代得到的近似值已经十分精确，所以可以舍去！

背后的故事

到这里，你已经基本上了解了卡马克快速平方根的算法原理，但其实，这个算法并不是卡马克所发明的！

卡马克当时担任的只是雷神之锤的 3D 引擎的开发者的一个角色，他只是代码的使用者，不过你在这个原始代码中看到的这段经典注释则确确实实出自于他手 —— what the fuck!

历史

1986年，一个名叫威廉卡汉（William Kahan）的人（后来也是IEEE754标准的制作者）用调整二进制数的方法写出来了求解平方根的算法，几年后，另一个程序员克莱佛莫勒（Cleve Moler）（MATLAB 背后的 MathWorks 公司创始人）发现了这段代码并介绍给了自己的同时格雷格沃什（Greg Walsh），而沃什在这段代码的基础上就写出了平方根倒数快速算法。

而沃什所在公司的另一位程序员加里塔罗里（Gary Tarol）看到了这段代码并创立了图形芯片公司 3dfx（90年代是可以碾压英伟达的存在）。塔罗里在 3dfx有一个同事布莱恩胡克（Brian Hook）跳槽了后来的雷神之锤（卡马克担任总负责人）的公司 id Software，便有了如今的卡马克平方根倒数快速算法。

一场比赛

在这段代码中不是有一个奇怪的值 0x5f3759df 嘛，普渡大学的数学家 Chris Lomont 看了以后觉得有趣，决定要研究一下卡马克的这个猜测值有什么奥秘。Lomont 也是个牛人，在精心研究之后从理论上也推导出一个最佳猜测值，和卡马克的数字非常接近, 0x5f37642f。Lomont 计算出结果以后非常满意，于是拿自己计算出的最佳猜测值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确的求得平方根。结果是卡马克还是赢。

最后 Lomont 怒了，采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数字要好上那么一丁点的数字，虽然实际上这两个数字所产生的结果非常近似，这个暴力得出的数字是 0x5f375a86。

如今

事实上，如果你现在使用 windows 等现代化操作系统以及现代化硬件来比较 sqrt 与快速平方根倒数的算法快慢，会发现两者相差甚微。

其中的原因是如今的底层硬件已经对 math.h 中的各种实现进行了底层优化，导致使用 sqrt 的算法速度也是非常快的。但这并不妨碍我们去了解历史上影响力最大的一段小代码片段的故事！这不就是计算机的乐趣嘛！: )