引言
在本文中，我们将深入研究快速计算平方根倒数，并了解神秘数字 0x5f3759df 的来源。 该算法在 id Software 开源 Quake III 引擎后出名。 在此过程中，我们还将学习浮点数和牛顿法。
Introduction
2005年，id 公司发布了《雷神之锤III：竞技场》的源代码。在源码中，游戏粉丝们找到了一个十分巧妙算法，很快出名了，该算法的唯一用途是计算平方根的倒数：
要我写计算平方根倒数的代码的话，我会这么写：
此处我用 C 语言，和《雷神之锤III》的编程语言相同。不过说真的，我没有自行实现开平方根的算法，那些对 C 语言底层和 CPU 设计更熟悉的人们已经解决了如何用电脑计算平方根的问题了，并且将该算法写进了 <math.h> 头文件，我们程序员直接将该头文件包含进程序就行了：
那么，怎么《雷神之锤III》的算法就很有趣了呢？当时 id 公司是怎么计算平方根倒数的呢？
乍看上去，完全解释不通啊。这个奇怪的 0x5f3759df 怎么冒出来的？算平方根的时候，这步有什么用途？为啥注释者在第二个注释那里极致享受了一下？在本文我将展示计算平方根时候使用到的一些很厉害的位操作，以及使用该方法的算法，名称叫平方根倒数速算法。
Why Care
首先，为啥游戏引擎要计算平方根的倒数呢？要在游戏引擎内实现物理效果、光影效果、反射效果的话要对使用向量进行单位化，因而我们需要计算平方根倒数。
要不然的话，向量会过短或过长。若在此基础上处理物理效果的话，会出错的。
众所周知，向量长度是 $x^2+y^2+z^2$ 的平方根。
如果你没反应过来，我保证你在二维空间里看到过类似的，就是勾股定理啊。
所以，如果你要将向量单位化（把向量长度缩放为1）的话，我们要将向量的三个维度都按照模长缩放长度。显然，如果我们将向量模长除以向量模长，结果显然为1：
那么，剩下的事情是将 x、y、z 轴的长度除以向量模长：
或者类似地，乘以向量模长的倒数：
你大概看明白了这其中的问题所在了吧，计算 $x^2+y^2+z^2$ 很简单，更重要的是真的很快。写代码时，你可以这么写：x*x + y*y + z*z：
这就是三个乘法相加啊。加法和乘法是计算机中的基本操作，而且被设计得十分快，然而平方根计算却是非常的慢，除法计算也快不到哪里去。
我们要处理上千的平面，每个向量都要单位化的话，这根本不行啊。但这也意味着，我们有机会可以提高该算法的速度，退一步讲，如果我们可以找到 x 的平方根倒数的近似结果，只要算法很快我们就能节省很多的时间。
平方根倒数速算法算的就是近似值，大约有 1% 的误差，不过速度是刚才那个的三倍：
The Code
再次查看代码，我们可以看到，代码开头基本没啥用，我们有一个形参 number 作为输入值，我们要计算这个数的平方根倒数。
首先有一个 32 位整型变量 i，我们还声明了两个 32 位的浮点数 x2 和 y。然后我们将 1.5 这个数存入一个常量 threehalfs：
下一行只是很简单地将形参的数值除以 2 后，存入 x2，完整的形参则放入 y。在此之后，魔法就发生了。
再看一下代码，好吧，你看这个代码的时间越长，这代码越解释不通，而代码右边的注释也帮不上什么忙。但是这至少提醒了我们，这个算法有三步，将这三步搞明白，算法就能向我们展示出它的奥妙了。
IEEE 754
在开始这三步之前，我们先来了解一下二进制数字。
代码第一行，我们声明了一个 32 位的整型变量，在 C 语言中，这被称为 long，也就是说我们有 32 个二进制位，我们可以利用这些位表示数字：
我觉得，你们肯定知道该如何表示吧。而在下一行我们声明了两个浮点数，在 C 中被称为 float，当然，我们也要用 32 个二进制位表示小数，你会怎么表示呢？
如果咱们一起设计小数的话，也许这是我们的设计方案之一：直接在中间放小数点就行了啊。在小数点之前，我们按照一般方法数出整数部分；而在小数部分之后，不出所料也是这样，不过这是二进制，所以我们不数十分位、百分位、千分位，我们数二分之一、四分之一、八分之一、十六分之一以及其所有的组合。
但这个想法实际上很糟糕，我们大幅缩小了这些位数可以表示的数字。之前我们可以表示大约 20 亿的数字，现在只有 3 万 2 千（$2^{15}=32768$）左右了。
幸亏，那些比我们聪明的人们想出个更好的办法来充分利用这 32 位，他们是从科学记数法中获得的灵感。这样，我们就能系统地表示数字了，诸如 $23000 = 2.310^4$，$0.0034 = 3.410^{-3}$：
在二进制体系中，我们也可以这么表示小数。举例，二进制数 11000 可以用 $1.1*2^4$ 表示：
规定此表示方法的标准称为 IEEE 754，和之前一样我们有 32 位，第一位是正负号位，如果该位是 0，则该数为正数；如果该位是 1，则该数为负数：
但是在游戏中，传给快速平方根倒数函数的值一定是正的。所以对于本文的剩余部分，我们不考虑正负号位，默认其为 0。
然后，接下来的 8 位决定指数部分，也就是 2 的 1 次方、2 次方、3 次方……以此类推：
我们可以利用这 8 位，表示介于 0 和 255 之间的数。但这个仍不是我们想要的，我们还想要负指数，这就是为什么这里出现的数用的时候得减去 127：
所以说我们实际得到的是 $2^{4-127}$ 而不是 $2^4$。如果我们需要的指数必须是 4 的话，这里必须要被设为 131，因为 131 - 127 = 4：
剩下的 23 位表示尾数。在科学记数法中，我们要先表示一位整数部分，然后是小数点，以及之后的小数部分。利用这 23 位，我们能表示从 0 到 $2^{23}$ 之间的任意一数：
当然，这还不是我们想要的。对于十进制科学记数法，我们需要尾数的表示范围在 1 到 10 之间，而在二进制记数法中则是 1 到 2 之间。
现在我们可以做一些之前就想做的事情了，在第一位的后面加小数点，这将会自动地使尾数的范围在 1 和 2 之间。但如此朴素的想法是对位数的浪费，你看啊，那些设计 754 标准的人们意识到在二进制里会出现一些在其他进制系统里不可能发生的事情。
看一下第一位，在科学记数法中，第一位按照定义一定不是 0，而在二进制中只有一个数是非零数：1。
要是第一位永远是 1 的话，我们就没有必要存储它了。这样我们节省了一位的存储空间，可以直接把小数点移动到左面，然后在回来转换的时候整数部分改成 1 即可。
现在尾数的表示区间是在 1 和 2 之间了。尽管 23 位能够表示 0 到 $2^{23}$ 之间的数，我们将表示数字的范围缩小到 0 和 1 之间，然后我们再补一个 1，就让其范围为 1 到 2 之间。
这就是 IEEE 754 标准的主要内容了，不过呢刚才说的都是“规格化数”，懂得多的人知道，这个标准还有非规格化数、非数、无穷和两个零（指 +0 和 -0）。
这里我们不会涉猎，因为在《雷神之锤III》，这些数恰巧是不会输入进算法的。
Bits and Numbers
对于本算法思考尾数、指数和二进制数字之间的关系十分有用。
我们有两个数，一个是尾数，一个是指数，分别用 23 个二进制位和 8 个二进制位表示。我们可以利用这些位数表示它们：
仔细想的话，是 $2^{23}$ 乘以指数 E 再加上尾数 M，因为将指数 E 乘以 $2^{23}$，相当于将其进 23 个二进制位：
所以，我们可以用这些二进制位表示浮点数了。我们可以用这个公式算出这些位数表示的实际数字，你应该很熟悉这个公式了：
这里，我们的指数被减去了 127，我们的尾数前面也多了一个 1。不过，现在的事情有点不同了，这里我们将完全没有缘由地对该公式取对数：
既然我们干的是计算机的事情，我们对数的底数取为 2。这里我们尽量简化，将指数取出，然后，我们卡住了，然而《雷神之锤》的开发者并没有。
开发者 Gary Tarolli 知道去掉对数的技巧（注：就现在所知，平方根倒数速算法最早由 Gary Tarolli 在 SGI Indigo 的开发中使用），技巧就是 $\log(1 + x)$ 的极限：
对于很小的 x，$\log(1+x)$ 大致相当于 x。仔细想的话，这个近似实际上是正确的：
但如果我们再加入一个校正系数称为 $\mu$，该校正系数可任意选择。若校正系数为 0，该近似函数将于 0 和 1 处与原函数所得因变量相同，不过最终得出结论，将校正系数 $\mu$ 设为该值，在 0 到 1 的范围内总体偏差最小：
回到公式，我们利用这个技巧。由于 M 已经被 $2^{23}$ 除去，它一定是 0 到 1 之间的数，再稍微整理一下，我们终于知道为啥我们要如此计算了：
这就是浮点数的二进制表示嘛。接下来，我们再想想我们之前干了啥：我们把公式取了对数，得到了二进制表示的形式，仅仅是缩放一下，移动了几位。在某种意义上，浮点数的二进制表示的就是它的对数形式。
懂了这么多，我们终于可以开始了解快速平方根倒数算法的三部曲了。
Evil Bit Hack
第一步实际上并不复杂，看起来复杂是因为这行代码是在内存地址方面做了些手脚。我们将数字存入 y，然后我们想做些很厉害的位操作。
然而浮点数没有做位操作的办法。为啥浮点数不能做位操作呢？因为它根本就不是为位操作设计的啊，浮点数本质上严格遵守 IEEE 754 标准，然而长整型数可以用来做位操作：
此处用一个小技巧举例，向左位平移能将这个数加倍，向右平移能将这个数减半：
当然，如果是个奇数的话最终会有舍去。不过呢，我们可以接受这点不准确，只要这个算法够快就好。
像其他的编程语言一样，C 语言提供了一种方法能将浮点数转换为长整型数，这种转换能满足程序员的大部分需要，就是将一个浮点数尽力转换为一个整型数。如果浮点数是 3.33 的话，它能将其转换为整型数也就是 3:
不过这个转换不是我们想要的。首先，我们对转换后的整型结果不感兴趣，我们想让其保持浮点数的形式；其次，浮点数背后的那堆二进制数也都被打乱了，我们不希望转换对二进制位有任何影响：
我们想要的是逐位将其转换到长整型形式。你要这么干的话，就要转换该数对应的内存地址，而不是数字本身：
首先我们得到了浮点数 y 的地址，这个是一个 float 类型的地址。然后我们将该地址转换成 long 长整型的地址：
地址本身没变，但是 C 语言认为该地址内的数是个长整型了。然后就读取该地址内的数字了，由于现在 C 语言认为该地址的类型是 long，它将按照长整型类型读取该地址里的数字。
像这样我们避开数字本身，利用该数所在的地址，运用了 C 语言的一些手段实现目的。这样该浮点数的二进制表示存入了 i，我们进入下一步。
WTF
算法的第二步有一点很显然，还记得吗，向左平移能将这个数加倍、向右平移能将这个数减半，但要对于指数进行该变换，会发生什么呢？
指数加倍，相当于对该数平方；指数减半，相当于对该数进行开方。现在对指数取负，就是将 x 开放后取倒数，这正是我们想要的。
所以我们要明白我们的目标就是这个，我们的数字存于 y，我们的目标就是计算 y 的平方根倒数。之前提到，直接计算相当困难，耗时很久，但是我们已经提取出 y 的二进制表示，我们也知道按照 IEEE 754 标准浮点数的二进制表示在某种意义上，就是该二进制的对数。
也就是说在 i 里我们存有 y 的对数，当然还有一些放缩和变换。我觉得要是处理对数的话，问题会变得简单。与其花很多力气直接计算 y 的平方根倒数，我们直接计算 y 的平方根倒数的对数，这里重新写成 y 的二分之一次方的对，这样指数可以提到前面，算这个可就简单多了。
你也许会想：“天哪，这里可有除法。刚开始不是说除法很慢吗？”。没错，不过这里我们可以进行位平移了，我们不算除以 2，直接将其向右平移一位，这就解释了为啥我们要算 i 的负数并将其向右移动一位。
不过， 0x5f3759df 这个数在这里有啥用啊？根据公式，我们还要放缩和位平移，我们来计算一下这个数从哪里来。
设 $\Gamma$ 为解，则解 $\Gamma$ 和 y 对数的负二分之一方相同，也就是 y 的对数乘以负二分之一。现在我们将该对数用二进制位来表示。
接下来我们仅需算出解 $\Gamma$ 的二进制位表示，此处略去细节，结果如下：
那个魔力数字，原来是按照校正系数 $\mu$ 放缩系数和平移系数来纠正误差的。现在我们有结果的位表示了，我们仅需再按照邪恶的位操作，将二进制位转换回浮点数结果即可。
实际上这不是确切结果，仅仅是一个近似值，这就是为啥我们需要第三步了。
Newton
上一步之后，我们得到了一个还可以的近似值了，但我们的结果或多或少有一些误差。按照牛顿迭代法，我们可以把还可以的近似值变成很好的近似值。
牛顿迭代法是一个用来寻找方程根的方法，也就是找到一个 x 使 f(x) = 0。它的实现方法是通过取近似的方式来得到一个更好的近似，通常你可以重复该过程，直到与实际结果十分接近：
然而实际上，我们已经很接近实际结果了。一次迭代得到的结果会有 1% 的误差。
牛顿迭代法仅需两件东西：方程对应函数和函数的导数。首先取一个 x 值，然后通过计算 f(x) 以及其导数来估计该值与实际根的偏差。我们将 f(x) 写为 y，导数则是 dy/dx，我们得到了 y 的偏差和 x 的偏差，以及 y 本身：
要得到 x 偏差，我们仅需将 y 除以其导数：
然后我们减去偏差，得到一个新的 x：
某些学霸现在意识到最后一行就是按照公式实现的牛顿迭代法，注意，作为平方根的 y 与 x 的平方根倒数是相同的。
我强烈建议大家自行验证该代码的最后一行，因为尽管该公式以及牛顿迭代法本身有除法，但代码里没有，这意味着该算法一直会很高速地运行。
Summary
现在我们终于搞明白快速平方根倒数算法了，仅仅是利用到了 IEEE 754 标准知识、C 语言特性的灵活利用、具有魔力的位平移以及牛顿迭代法背后的运算规则。