引擎技术漫游指南:卡马克的魔法


3楼猫 发布时间:2023-10-26 23:32:28 作者:TirelessKoriel Language

我们先从一个简单的问题入手。
在二维平面中有一个向量,它代表了我们生活在二维空间中的主角当下运动的状态。向量的方向即运动方向,向量的大小即速度大小。如果我们不在于主角的移动速度,只是想得到他的运动方向,也就是我们常说的方向向量,我们该如何计算?
很简单,我们利用勾股定理 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 - 127M / 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) 的专用硬件顶点着色器中实现。

© 2022 3楼猫 下载APP 站点地图 广告合作:asmrly666@gmail.com