本文系用户投稿，不代表机核网观点
⚠️ 未经作者授权 禁止转载
我们先从一个简单的问题入手。
在二维平面中有一个向量，它代表了我们生活在二维空间中的主角当下运动的状态。向量的方向即运动方向，向量的大小即速度大小。如果我们不在于主角的移动速度，只是想得到他的运动方向，也就是我们常说的方向向量，我们该如何计算？
很简单，我们利用勾股定理 a^2 + b^2 = c^2 即可求得向量长度的二次方，随后我们再对结果做平方根运算，得到的值为向量的大小，最后我们再用原向量除以向量大小，得到的结果即为方向向量，这个计算方法通常被叫做归一化（normalize），得到的向量我们也可以叫做单位向量。
方向向量不仅是在游戏中的物理模拟非常重要，对于渲染来说也是一样的。它能用来指明光线与三角面的方向，而方向则决定了光照程度、材质质感等关键效果的计算结果。
我们用简单的式子来描述一下这个计算过程：V / ||V||，其中 V 为任意向量，||V|| 为向量大小。也就是说，我们只要得到一个向量大小的倒数，即 1 / ||V||，就可以快速计算出方向向量。而 ||V|| 又涉及到平方根操作，所以最终我们要计算的其实是平方根的倒数。
当然我们可以把倒数与平方根运算分开操作，但是在游戏开发的蛮荒时代，平方根运算并不容易，倒数的操作也不是非常快速的简单运算。被机能困在旧时代的开发者们渴望得到一个快速算法，能够解决游戏开发中涉及到的大量的平方根倒数运算。
那么这个问题该怎么解决呢？
牛顿的咒语
牛顿的咒语
平方根该怎么计算？当然你不能使用计算器。
举个例子来说，我们就求 √2 吧。
我们可以做如下推导：
x = √2
=> x^2 = 2
=> 0 = x^2 - 2
所以我们可知，函数 f(x) = x^2 - 2 与 x 轴的交点即为 √2 的值，函数图像如下：
观察图像，我们可以推测这个值介于 1 与 1.5 之间。
上面的函数与方程只能让我们得到 √2 这种形式的数值表达形式，那么我们怎么求得 √2 具体值呢？
如上图所示，我们画一条线，这条线与函数 f(x) 的曲线相切于点（2,2），并且这条线与 x 轴的交点十分的接近我们要求得的那个点，这是一个好现象。也就是说，我们只要找出这条线的函数形式，并求得这条线与 x 轴的交点，我们就能的到一个 √2 的近似值。
那么我们怎么求得这条线的函数呢？
首先我们知道这条线经过点（2,2），那么如果我们还能得到线的斜率，通过点斜式即可求得函数。
因为这条线本质上是函数 f(x) 在点（2,2）切线，所以它的斜率就是 f(x) 在这一点的导数。
所谓的函数的导数，它表示这个函数在某一点上变化的趋势与程度。这一概念其实与经过切点的切线是保持一致的。切线的倾斜程度（即斜率）不就是代表了这条线上的值的变化程度吗？
我们用 f '(x) 来表示函数 f(x) 的导函数，对于 f(x) = x^2 - 2 来说，f '(x) = 2x。
其实对于类似 f(x) = ax^n + b 这种形式的函数来说，导函数的形式非常简单，就是 x 的指数减一，得到的值与 x 的系数相乘，偏移量 b 因为是常数，所以会被省略。按照这个规则，f(x) = ax^n + b 的导数就是 f ‘(x) = an * x^(n-1)，所以 f(x) = x^2 - 2 导数即为 f '(x) = 2x。
我们将点（2,2）带入函数 f '(x) = 2x 求得此处的导数为 4，即经过此处的切线的斜率为 4。
现在切线的斜率与切线经过的点都有了，带入点斜式我们可得：
y - b = k(x - a)；a，b 为函数图像经过的点
y - 2 = 4(x - 2);
y = 4x - 6;
即切线函数为 g(x) = 4x - 6，并计算 0 = 4x - 6，所以切线与 x 轴的交点为（1.5, 0）。那么我们要求的 √2 近似值即为1.5。
1.5 这个值明显并不精确，它有点过于“近似”了。从图上我们也能观察到它与真正的 √2 的值还是有一定差距的。但是我们已经掌握了求近似值的方法，那么我们是否可以重复上述过程，来进一步接近正确值呢？
答案当然是肯定的。我们已经得到了一个接近 √2 的值 1.5，那我们可以尝试以这个值为基础，再做一次上面的运算。
首先我们将 1.5 带入 f(x)，求得一个新的切点：
f(x) = x^2 - 2 = 1.5^2 - 2 = 0.25
新切点即为（1.5, 0.25）。
再导入导函数求切线斜率：
f '(x) = 2x = 2 * 1.5 = 3
最后带入点斜式，得到切线函数：
y - b = k(x - a)；
y - 0.25 = 3(x - 1.5);
y = 3x - 4.25;
3x = 4.25;
x ≈ 1.41666
我们也可以从图像中观察到这个结果：
我们将这几条线与 x 轴相交的部分放大观察：
正如我们设想的那样，运算的得到的结果作为下一次运算的初始值，递归运算下去，我们就能不断的接近那个我们想要的值。
我们再从全局的角度观察一下我们的计算的流程。我们有两个关键函数，f(x) 和它的导函数 f '(x)，用它们我们可以求得切点与斜率，再带入点斜式求得切线函数，最后得到切线与 x 轴的交点，即为我们想要的值。从这个过程来看，我们好像可以总结出一个更简单的公式，来表示这个计算过程，我们来尝试一下。
点斜式 y - b = k(x - a)，因为我们最后要求得的值为函数与 x 轴交点中的横坐标值，即点斜式中的 a 值；又因为 x 轴交点的纵坐标必为 0，所以我们将点斜式变化为：
a = x - y/k
其中 x 与 y 组成我们当前已知的切点（例如上述第一次计算过程中的点 2,2），y 由 f(x) 计算得到，所以 y = f(x)；k 为斜率，由导函数 f ’(x) 计算的到，所以 k = f ’(x)。
我们还知道 a 为我们要求得的最新的近似值，我们设它为Xn+1；x 为我们当前已经求得的近似值，我们设它为Xn。
将上述各个值的转换导入式子，我们可得：
Xn+1 = Xn - f(Xn) / f '(Xn)；
带入具体的实例测试一下？
第一个切点我们选（2,2），并且我们有 f(x) = x^2 - 2，f '(x) = 2x，带入公式即为：
Xn+1 = 2 - (2^2 - 2)/(2*2) = 2 - 1/2 = 1.5。
恭喜！我们又一次发明了牛顿的咒语：牛顿迭代法。
在这个过程中，我们不难发现牛顿灵活的运用了近似值这个概念，以近似值的方式，不断逼近精确值。与其说是算法， 不如说是一种策略。
对于这条所谓的咒语来说，它不仅简单，而且还特别适合用代码实现，我们可以用 for 或 while 为循环计算提供动力，直到计算出我们想要的精度，一次平方根计算就完成了。
但是，我们开头也提到过，这是一个蛮荒的时代，计算机的算力有限，而我们要完成的却有很多，我们很难使用一个依赖循环次数得到确切值的算法。循环就像一个吞噬算力的猛兽，它很快就会让机器的计算能力变得捉襟见肘。
如果只需要念一次咒语就好了。
没错，只进行一次牛顿迭代法。不是不可能，只不过需要我们首次提供的 Xn 就足够精确罢了，足够精确的 Xn 得到更精确的 Xn+1。
这好像是个死循环。
如果卡马克的魔法没有诞生的话。
卡马克的魔法
卡马克的魔法
假设我们有一个数字，30000，用科学计数法表示它即为，3 * 10^4。不难发现 10 的指数其实就是数值中 0 的个数。
十进制的规则深入我们的骨髓，其中的道理也无需多言。
那么二进制呢？
假设我们有一个二进制数 10000（十进制值为 16），我们类比上面的十进制规则，那么二进制 10000 的科学计数法表示形式即为 1 * 2^4。此处与十进制不同的是，二进制是以 2 为底的指数形式。你可以仔细品味一下乘法的本质，再从十进制类比到二进制。
这一切都说得通，而且非常自然。
我们再进一步。
十进制数 31000，用科学计数法表示即为 3.1 * 10^4，观察这个表达形式，我们不难发现，10 的指数本质其实是在表明小数点的移动位数。
那么类比到二进制中，二进制数 11000 （十进制值为 24）的科学计数法表达形式即为 1.1 * 2^4，符合我们上面观察到的规律。
这就是二进制的位移规则。
如上图所示，浮点数在计算机中的二进制存储模式被分为三个部分，首位是符号位，用来标识数字的正负号；后八位是偏执指数位，用来表示在科学计数法下小数点偏移的位数；最后的二十三为则表示在科学计数法下，小数点后面的数字排列，这部分融合了原数字的整数与小数两个部分。
我们以一个具体的数字来举例，二进制 111 的科学计数法表达形式为：
1.11 * 2^2
以浮点数编码规则编码后，在计算机中的二进制形式为：
其中指数位的计算规则是，小数点偏移量 + 127。以 111 为例，科学计数法表达后小数点向左偏移两位，所以指数位为 2 + 127 = 129，二进制表达即为 10000001。
至于尾数位，因为二进制科学计数法的规则就是小数点会停留在首个不为 0 的位数后面，所以我们只需要去掉科学计数法中的“1.”，剩下的部分就是尾数位的本体。
以上就是我们对浮点数编码规则的理解。
关于浮点数，我们前面所讨论的一直都是科学计数法到浮点数二进制形式的转换，那我们就忍不住好奇，将这个过程反过来呢？浮点数二进制形式如何转化为科学计数法呢？
很简单，对于偏置指数位来说，我们只需要用偏执指数位 - 127，结果就是科学计数法中的小数点偏移量，就上述中的例子来说：
100000001 - 01111111
= 129 - 127
= 2
对于尾数位来说，反向推到科学计数法的表达形式时，需要在尾数位前面加上“1.”，就上述中的例子来说，也就是在 “11000……” 前面加上 “1.”，即 1.11，0 被省略。
这还不够，根据上面的文字描述，我们可以试着总结出一条公式。
我们设尾数位为 M，因为尾数位共 23 位，尾数位反推至科学计数法时，小数点要出现在尾数位形式之前，那么根据我们前面提到的二进制位移规则，尾数位小数点向左移动 23 位，即 M / 2^23，最后别忘了，我们还要加上一个“1.”，所以最终公式应该为：1 + M / 2^23；
接着我们设偏置指数位为E，则科学计数法中小数点偏移量的计算为：E - 127。
结合尾数位部分的公式，我们可以得到（这里我们省略的符号位，因为加入符号位很简单，而且后续我们的讨论也不需要符号位）：
科学计数法表达形式 = (1 + M / 2^23) * 2^(E - 127)
对比一下例子中的科学计数法表达：1.11 * 2^2
那么这个公式有什么用呢？我们先来随意摆弄一下它看看会发生什么。
因为公式很明显被乘法分为两个部分，并且其中一部分还是 2^(E - 127) 这种 2 的指数形式，那么我们不防试试对数？
对数的运算法则很有意思：
它能将乘法的两部分，变为加法（第一条）；还能将指数转变为乘法（第四条）。这两条都存在于我们的公式之中，那我们就来试试。
设 a = (1 + M / 2^23) * 2^(E - 127)，我们对两边都进行以 2 为底的对数运算，可得：
㏒₂a = ㏒₂((1 + M / 2^23) * 2^(E - 127))
㏒₂a = ㏒₂(1 + M / 2^23) + ㏒₂2^(E - 127)
㏒₂a = ㏒₂(1 + M / 2^23) + ㏒₂2 * E - 127
㏒₂a = ㏒₂(1 + M / 2^23) + E - 127
唉？有趣，2 的指数形式消失了，只剩下 E - 127，从某种角度来说，整体运算变得简单了。那 ㏒₂(1 + M / 2^23) 部分呢？
我们知道 M 作为尾数位的二进制形式，具体二进制取决于我们浮点数的具体数值，所以这里可以把它看做是一个变量，M / 2^23 自然也是随着 M 的变化而变化的，我们不防暂时先将他们看做一个整体 x。于是 ㏒₂(1 + M / 2^23) 就变成了 ㏒₂(1 + x)。这就又简化了我们的式子，因为 ㏒₂x 本身就是一个很常见的对数函数，图像如下：
而 ㏒₂(1 + x) 只需对图中图像稍作位移：
这还没完，因为我们知道 M / 2^23 必然大于 0 小于 1，所以 ㏒₂(1 + x) 中，x 的取值范围也在 0 至 1 之间，那么我们放大这个图像的细节到 x 的 0 至 1 区间：
从图像中可以看到，这条曲线在 0 到 1 之间的部分，其实近乎于直线，再回想一下我们从牛顿的咒语中学到了什么？
没错，近似值。在我们实际的工程应用中，很多数值都不是十分精确的，或者说对数值的精确度要求是变化的。所以既然它接近于一条支线，那我们为什么不就把它当做一条支线呢？
因为曲线经过点原点与点（1，1），所以我们可以把它看做 y = x：
从图上可以观察到，y = x 与 ㏒₂(1 + x) 在 0 至 1 区间的值，真的非常相近。
那么，因为：
㏒₂(1 + M / 2^23) ≈ M / 2^23
所以：
㏒₂a = ㏒₂(1 + M / 2^23) + E - 127
≈ M / 2^23 + E - 127
到这里还不够，我们再稍作变形，将 M / 2^23 + E - 127 中 M / 2^23 + E 提取 1 / 2^23 可得：
1 / 2^23 * (M + E * 2^23) - 127
最后我们终于来到了这里。
观察式子中的 M + E * 2^23，你看它像什么？
还是 1.11 * 2^2 的二进制形式为例，它的二进制形式为：
其中，M = 11000000000000000000000；
E = 10000001；
那么，E * 2^23 根据二进制位移规则，就是将小数点向右移动 23 位，就是说要在 E 后面添加 23 个 0：
E * 2^23 = 10000001 00000000000000000000000；
于是可得：
M + E * 2^23
=
11000000000000000000000
+
10000001 00000000000000000000000
=
10000001 11000000000000000000000
原来，M + E * 2^23 就是这个 M 和 E 所代表的那个数字的二进制形式。
那么根据我们的推导过程：
㏒₂a = ㏒₂(1 + M / 2^23) + E - 127
≈ M / 2^23 + E - 127
= 1 / 2^23 * (M + E * 2^23) - 127
我们可以得到一个结论：
对于数字 a，设 A 为浮点数 a 二进制形式，那么：
㏒₂a ≈ 1 / 2^23 * A - 127
也就是说我们得到了一个快速计算对数的算法，只要利用浮点数二进制就可以快速计算。
还记得我们最开始的目的吗？我们是为了计算平方根倒数，那这个快速计算对数的方法和平方根倒数计算又有什么关系呢？别着急，我们接着来。
我们设 a 的平方根倒数，即 1/√a ，为 y，那么我们有：
㏒₂y = ㏒₂a^(-1/2)
= -1/2 * ㏒₂a = -1/2 * ((1 / 2^23) * A - 127)
也就是说：
㏒₂y = -1/2 * ((1 / 2^23) * A - 127)
y 虽然是 1/√a 的符号表现，但是，y 最终一定有一个值对不对？它既然是一个值，那就必然可以表达为浮点数的二进制形式，我们就设浮点数 y 的二进制形式为 Y。这个时候我们就可以 ㏒₂y 也转化为二进制快速计算对数的形式：
㏒₂y = (1 / 2^23) * Y - 127
而 Y 是浮点数 y 的二进制形式，也就代表了 y 的值，即 1/√a，也就是我们想要求得的 a 的平方根倒数。
终于，一切都说的通了：
(1 / 2^23) * Y - 127 = -1/2 * ((1 / 2^23) * A - 127)
Y = 381 * 2^22 - 1/2 * A
Y = 1598029824 - 1/2 * A
如果我们用符号 “>>” 表示二进制向右的移动，那么考虑到我们前面所提到的二进制位移规则， 1/2 * A 就等同于 A >> 1，它表示将二进制 A 向右移动一位；如果我们将 1598029824 转化为十六进制，那么最终的式子就变为：
Y = 5F400000 - A >> 1
根据最终我们得到得的这个式子，我们得出一个结论：
将浮点数 a 的二进制形式所代表的值乘以0.5，再与数值 1598029824 作差，结果的绝对值即为 A 的平方根倒数。
这就是卡马克在雷神之锤 3 中所写下的快速平方根倒数算法。它巧妙的将浮点数的编码规则与平方根倒数的计算相结合，一切都显得如此自然又仿佛是巧合。
我们看一下具体的代码实现：
其中第七行代码，先将 y 值转化为 long 型，因为整型的编码就是简单的二进制与十进制的转换，不涉及类似于浮点数编码规则对二进制数的特别解读，所以转化为 long 型，就相当于是原封不动的保存了浮点数 y 的二进制。
第八行代码与我们上面推导出的 Y = 5F400000 - A >> 1 相似，不同的点在于减数，代码中的减数 0x5f3759df 看起来比我们的 5F400000 要更加精确。这也是事实，还记得我们上面提到过的 y = x 与 ㏒₂(1 + x) 在 0 至 1 区间的图像对比吗？
我们能看到，虽然两个函数的图像在 0 至 1 区间内相似，但是归根结底还是不同的，然而我们的所有推导，都是建立在他们非常相似的基础上，它们越相似，我们的误差就越小，算法的精度就越高。如果我们将 y = x 变为 y = x + 0.05，则图像如下：
这个 0.05 我们可以叫做修正值。这个修正值在我们的计算过程中，最终会以常数的形式与 5F400000 相加，实际代码中的 0x5f3759df 就是这么来的。这个修正值可以按照一定的策略来做尝试，比如我们可以使用二分法，快速缩减尝试范围，最红确定一个相对精确的修正值。
最后第十一行代码，使用了牛顿迭代法进一精确我们得到的数值。具体代码的逻辑为：
Y = 1/√x
=> Y^2 = 1/x
=> 1/Y^2 - x = 0
即，我们有函数 f(Y) = 1/Y^2 - x，此函数当 x 为 0 时的图像如下：
由图可知，当 x = 0 时，函数图像无限接近 x 轴，但并不与 x 轴相交；如果 x >= 0，我们可以想象，函数图像会向下位移，之后与 x 轴相交处，即是我们要求得的 x 的平方根倒数。
因为 x 是我们想要求得平方根倒数的目标数值，所以对于每次计算来说，这都是一个常数，那么根据牛顿迭代法可得：
已知，f '(Y) = -2/Y^3
所以，
Yn+1 = Yn - f(Yn) / f '(Yn)
=> Yn - (1/Yn^2 - x) / (-2/Y^3)
=> Yn * (3/2 - x/2 * Yn^2)
此式子就对应代码中：
y = y * ( threeHalfs - ( x2 * y * y ) );
这与我们在牛顿的魔法一节最后提到的，以一个相对精确值得到更精确值的目标相对应。
最后值得一说的是，其实这个快速计算平方根倒数的算法并非出自卡马克之手，不然卡马克的注释也不会是 WTF？了。
对于这个算法最初的发明者，我们可以参考 wiki 上的说法：
威廉·卡汉 (William Kahan) 和 K.C. 伯克利分校的 Ng 于 1986 年 5 月写了一篇未发表的论文，描述了如何通过使用比特微调技术（bit-fiddling techniques）和牛顿迭代法来计算平方根。 20 世纪 80 年代末，Ardent Computer 的 Cleve Moler 了解了这项技术，并将其传授给了他的同事 Greg Walsh。 Greg Walsh 设计了现在著名的常数和快速反平方根算法。 Gary Tarolli 正在为 Kubota（该公司当时为 Ardent 提供资金）提供咨询，并可能在 1994 年左右将该算法引入 3dfx Interactive Jim Blinn 在 1997 年 IEEE 计算机图形和应用专栏中演示了平方根倒数的简单近似值计算。对其他当代 3D 视频游戏的逆向工程中，发现了动视 1997 年 Interstate '76 中算法的变体。 Quake III Arena 是一款第一人称射击游戏，由 id Software 于 1999 年发布，并使用了该算法。 Brian Hook 可能将 3dfx 的算法引入了 id Software。2000 年，中国开发者论坛 CSDN 上出现了对该代码的讨论，Usenet 和 gamedev.net 论坛于 2002 年和 2003 年广泛传播了该代码。人们开始猜测谁编写了该算法以及该常数是如何导出的。 有些人猜测是约翰·卡马克。Quake III 的完整源代码在 QuakeCon 2005 上发布，但没有提供答案。 2006 年，当 Greg Walsh 联系 Beyond3D 时，作者身份问题得到了解答，因为他们的猜测在 Slashdot 上广受欢迎。 2007 年，该算法在一些使用现场可编程门阵列 (FPGA) 的专用硬件顶点着色器中实现。
评论区
共 4 条评论热门最新