《雷神之锤III竞技场》与快速逆平方根算法

作者： Paul King

发布时间：2023-02-28 上午12:05

简介

1999 年，id Software 发布了多人第一人称射击游戏《雷神之锤III竞技场》。

Quake III

当游戏源代码向世界发布时，其中包含一个此前未知的算法，称为快速逆平方根。我们不打算深入解释该算法，但其当时的意义在于，它为方程 \$f(x) = 1/sqrt(x)\$ 提供了一个非常好的近似值，该方程用于向量归一化和游戏背后的其他数学运算，广泛用于包括 3D 着色在内的众多图形方面。与使用传统数学库相比，快速算法计算结果的速度快 3-4 倍，并且结果误差在 0.2% 以内。

这是代码：

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 f☼⁕k?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

它为什么看起来很奇怪？其中有许多部分是你不会期望出现在平方根计算中的。有一些技巧用于浮点数和长整型 IEEE 754 32 位表示之间的相互转换，一个“魔术”常数，一些位移，以及一些牛顿法的点缀，以确保良好的效果。

这些细节已在众多博客中详细解释，例如快速平方根近似和理解《雷神之锤》的快速逆平方根，以及视频，例如快速逆平方根 — 《雷神之锤III》算法和《雷神之锤III竞技场》的未知算法解释。甚至还有一些网站展示了多种语言（包括 Groovy）的算法。

游戏发布后不久，英特尔在 x86 处理器中添加了 SQRTSS 和 RSQRTSS 指令。以下是 JDK 增强功能之一，为 float Math.sqrt() 提供 SQRTSS 支持。后来的博客解释说，由于这些及其他硬件的进步，该算法现在主要只在历史传说层面具有相关性，例如《雷神之锤3》逆平方根技巧的消亡和快速逆平方根仍然快速吗？。

让我们看看几种计算逆平方根的替代方法，包括快速逆平方根算法的两种变体。我们将在 Java 和 Groovy 中进行此操作，并查看 `float` 和 `double` 实现。我们可以进一步尝试不同次数的牛顿法校正迭代，但我们将此作为读者的练习。😊

Java Float 实现

有许多地方可以看到该算法的 Java 示例，也有许多方法可以运行微基准测试。我们使用的代码与此 gist 略有相似。我们可以做得更精细，使用 jmh，但我们这里不打算进行全面的性能分析，我们只是想要一些粗略的估计数字。也许那是后续博客的主题。

与所有微基准测试一样，在对结果赋予过多的权重之前，您应该格外小心。在不同的机器上，使用不同的 Java 和 Groovy 版本等，数字将有所不同。我们使用了 Groovy 4.0.8 和 Java 17.0.2。

我们从最明显的实现开始，即使用 `Math.sqrt` 方法，该方法接受并返回 `double` 类型的值。

public static float invsqrt0(float x) {                   // Java
    return 1.0f / (float)Math.sqrt(x);
}

与上述 C 语言实现相比，Java 没有位类型转换技巧，但它有自己的方法来获取位表示。

public static float invsqrt1(float x) {                   // Java
    float xhalf = 0.5f * x;
    int i = Float.floatToIntBits(x);
    i = 0x5f3759df - (i >> 1);
    x = Float.intBitsToFloat(i);
    x = x * (1.5f - (xhalf * x * x));
    return x;
}

作为替代，我们可以使用字节缓冲区来回转换位。

private static ByteBuffer buf = ByteBuffer.allocateDirect(4);      // Java

public static float invsqrt2(float x) {
    float xhalf = 0.5f * x;
    int i = buf.putFloat(0, x).getInt(0);
    i = 0x5f3759df - (i >> 1);
    x = buf.putInt(0, i).getFloat(0);
    x = x * (1.5f - (xhalf * x * x));
    return x;
}

Java Double 实现

同样，我们将从最明显的实现开始。

public static double invsqrt0(double x) {                   // Java
    return 1.0d / Math.sqrt(x);
}

同样，使用内置方法获取位表示。

public static double invsqrt1(double x) {                   // Java
    double xhalf = 0.5d * x;
    long i = Double.doubleToLongBits(x);
    i = 0x5fe6ec85e7de30daL - (i >> 1);
    x = Double.longBitsToDouble(i);
    x *= (1.5d - xhalf * x * x);
    return x;
}

代码类似于 float 版本，但“魔术”常数已针对双精度浮点数更改。

字节缓冲区替代方案

private static ByteBuffer buf = ByteBuffer.allocateDirect(8);      // Java

public static double invsqrt2(double x) {
    double xhalf = 0.5d * x;
    long i = buf.putDouble(0, x).getLong(0);
    //long i = Double.doubleToLongBits(x);
    i = 0x5fe6ec85e7de30daL - (i >> 1);
    x = buf.putLong(0, i).getDouble(0);
    x *= (1.5d - xhalf * x * x);
    return x;
}

我们也可以比较一下 `Math.pow` 方法

public static double invsqrt3(double x) {                   // Java
    return Math.pow(x, -0.5d);
}

（我们也可以对 `float` 执行此操作，但这并不能为我们的分析增加太多内容，因为它无论如何都会调用此双精度方法。）

Groovy Float 实现

所有这些示例都启用了静态编译。我们需要速度，并且不进行任何元编程，所以不需要 Groovy 的动态能力。

我们的代码在显而易见的情况下与 Java 类似

static float invsqrt0(float x) {
    1.0f / Math.sqrt(x)
}

同样，对于快速算法，代码与 Java 类似。

static float invsqrt1(float x) {
    float xhalf = 0.5f * x
    int i = Float.floatToIntBits(x)
    i = 0x5f3759df - (i >> 1)
    x = Float.intBitsToFloat(i)
    x *= 1.5f - (xhalf * x * x)
}

再次使用字节缓冲区

private static ByteBuffer buf = ByteBuffer.allocateDirect(8)

static float invsqrt2(float x) {
    float xhalf = 0.5f * x
    int i = buf.putDouble(0, x).getInt(0)
    i = 0x5f3759df - (i >> 1)
    x = buf.putInt(0, i).getDouble(0)
    x *= 1.5f - (xhalf * x * x)
}

我们也可以尝试 Groovy 的 `**` 运算符（`power` 方法）

static float invsqrt4(float x) {
    (x ** -0.5).floatValue()
}

Groovy Double 实现

标准方法现在应该很熟悉了

static double invsqrt0(double x) {
    1.0d / Math.sqrt(x)
}

快速算法

static double invsqrt1(double x) {
    double xhalf = 0.5d * x
    long i = Double.doubleToLongBits(x)
    i = 0x5fe6ec85e7de30daL - (i >> 1)
    x = Double.longBitsToDouble(i)
    x *= (1.5d - xhalf * x * x)
}

结合字节缓冲区

private static ByteBuffer buf = ByteBuffer.allocateDirect(8)

static double invsqrt2(double x) {
    double xhalf = 0.5d * x
    long i = buf.putDouble(0, x).getLong(0)
    i = 0x5fe6ec85e7de30daL - (i >> 1)
    x = buf.putLong(0, i).getDouble(0)
    x *= (1.5d - xhalf * x * x)
}

使用 `Math.pow`

static double invsqrt3(double x) {
    Math.pow(x, -0.5d)
}

再次使用 Groovy 的 `**` 运算符（`power` 方法）

static double invsqrt4(double x) {
    (x ** -0.5).doubleValue()
}

结果

以下是执行上述方法的结果。我们使用了与之前提到的 gist 类似的一个测试框架，但我们计算了 100,000 个随机数的逆平方根，而不是 10,000 个，并且我们在计时循环中使用了 1000 次迭代，并计算了每 100,000 次计算的平均执行时间。

算法与实现（时间毫秒/秒）	Java 浮点数	Java 双精度浮点数	Groovy 浮点数	Groovy 双精度浮点数
Math.sqrt	0.216	0.254	0.359	0.245
快速逆平方根	0.230	0.236	0.357	0.127
带字节缓冲区的快速逆平方根	0.337	0.364	0.486	0.187
Math.pow(x, -0.5d)		8.949		8.997
x ** -0.5			0.737	1.807

分析

对于所有示例，使用字节缓冲区总是比原始算法慢，因此我们可以将其排除作为优化。我们也可以排除使用 `Math.pow` 方法。它比 `Math.sqrt` 慢得多。有趣的是，Groovy 的 `**` 运算符（`power` 方法）虽然仍然是较慢的选择，但明显优于 JDK 库的 `pow` 实现。

Java 快速算法对于浮点数较慢，对于双精度浮点数仅略微快一点。在大多数情况下，额外的代码复杂性带来的微小性能提升似乎不值得。

Groovy 的浮点实现稍微慢一些。这是因为 Groovy 在大多数计算中使用双精度浮点数，并在步骤之间将其转换回浮点数。这可能是未来优化的一个领域。

Groovy 的双精度实现至少与 Java 一样快。有趣的是，快速算法在 Groovy 中似乎更快。如果你确实需要速度，这些算法似乎值得进一步研究，但考虑到双精度浮点数的额外精度，你可能需要运行额外的牛顿法迭代，而这些迭代可能会抵消任何节省的时间。

有趣的是，对于所有实现，随机数快速算法的误差都非常接近 6E-5。

结论

我们使用 Java 和 Groovy 对快速逆平方根算法进行了微基准测试。结果呢？你可能再也不需要担心快速逆平方根算法了！但如果你真的有需求，它仍然相对较快，但你应该对你的应用程序进行基准测试，看看它是否真的有帮助。

《雷神之锤III竞技场》与快速逆平方根算法

简介

Java Float 实现

Java Double 实现

Groovy Float 实现

Groovy Double 实现

结果

分析

结论

参考文献