Quake III Arena 和快速平方根倒数算法

作者: Paul King
发布时间: 2023-02-28 12:05AM

介绍

1999 年，id Software 发布了 Quake III Arena，一款多人第一人称射击游戏。

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 位表示之间进行转换，一个“神奇”常数，一些位移位，以及一点牛顿法作为补充。

许多博客已经详细解释了这些细节，例如快速平方根近似和理解 Quake 的快速平方根倒数，以及视频，例如快速平方根倒数 - Quake III 算法和 Quake III Arena 的未知算法解释。甚至还有一些网站展示了该算法在多种语言中的实现（包括 Groovy）。

在游戏发布后不久，英特尔就在 x86 处理器中添加了 SQRTSS 和 RSQRTSS 指令。以下是一项 JDK 改进，为 float Math.sqrt() 提供 SQRTSS 支持。后来的博客解释说，由于这些指令和其他硬件进步，该算法现在主要只具有历史意义，例如 Quake 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 个，并且在计时循环中使用了 1,000 次迭代，并找到了每 100,000 次计算的平均执行时间。

算法 vs 实现（时间 m/s）	Java Float	Java Double	Groovy Float	Groovy Double
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 快速算法对于 float 来说更慢，而对于 double 来说只有轻微的提升。在大多数情况下，承担额外的代码复杂度似乎不值得获得微小的性能提升。

Groovy float 实现略微慢一些。这是因为 Groovy 使用 double 进行大多数计算，并在步骤之间转换回 float。这是未来可能优化的一个领域。

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

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

结论

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

Quake III Arena 和快速平方根倒数算法

介绍

Java Float 实现

Java Double 实现

Groovy Float 实现

Groovy Double 实现

结果

分析

结论

参考文献