在Rust中编写Quake的快速InvSqrt()函数是否可行？

Question

这只是为了满足我的好奇心。

这个有没有实现：

float InvSqrt (float x)
{
   float xhalf = 0.5f*x;
   int i = *(int*)&x;
   i = 0x5f3759df - (i>>1);
   x = *(float*)&i;
   x = x*(1.5f - xhalf*x*x);
   return x;
}

在 Rust 中有吗？如果有，请发代码。

我尝试了但失败了。我不知道如何使用整数格式编码浮点数。这是我尝试的代码：

fn main() {
    println!("Hello, world!");
    println!("sqrt1: {}, ",sqrt2(100f64));
}

fn sqrt1(x: f64) -> f64 {
    x.sqrt()
}

fn sqrt2(x: f64) -> f64 {
    let mut x = x;
    let xhalf = 0.5*x;
    let mut i = x as i64;
    println!("sqrt1: {}, ", i);

    i = 0x5f375a86 as i64 - (i>>1);

    x = i as f64;
    x = x*(1.5f64 - xhalf*x*x);
    1.0/x
}

参考资料：
1. Quake3的Fast InvSqrt()起源 - 第1页
2. 理解Quake的快速倒数平方根
3. 快速倒数平方根.pdf
4. 源代码：q_math.c#L552-L572

- Flyq

---

C#版本：是否有可能在C#中编写Quake的快速InvSqrt（）函数？ - Flyq

据我理解，由于违反了C严格别名规则，此代码在C语言中属于未定义行为。进行这种类型的别名操作的标准方法是使用“联合体”(union)。 - trent

我不认为union也行。memcpy绝对可以，尽管它很冗长。 - Matthieu M.

使用联合体进行类型转换是C语言中完全有效的，但在C++中是无效的。 - Salem

我认为从单纯好奇的角度来看，这个问题是可以的，但请了解时代已经改变了。在x86架构中，Pentium III于1999年推出的rsqrtss和rsqrtps指令比这段代码更快、更准确。ARM NEON有一个类似的指令叫做vrsqrte。而Quake III所使用的任何计算，现在可能都会在GPU上执行。 - benrg

Answer 1

我不知道如何使用整数格式对浮点数进行编码。

有一个函数可以做到这一点: f32::to_bits，它返回一个u32。还有另一个方向的函数:f32::from_bits，它以u32作为参数。这些函数比mem::transmute更可取，因为后者是不安全的，难以使用。

有了这个，下面是InvSqrt的实现:

fn inv_sqrt(x: f32) -> f32 {
    let i = x.to_bits();
    let i = 0x5f3759df - (i >> 1);
    let y = f32::from_bits(i);

    y * (1.5 - 0.5 * x * y * y)
}

(游乐场)

---

这个函数在x86-64上编译成以下汇编代码：

.LCPI0_0:
        .long   3204448256        ; f32 -0.5
.LCPI0_1:
        .long   1069547520        ; f32  1.5
example::inv_sqrt:
        movd    eax, xmm0
        shr     eax                   ; i << 1
        mov     ecx, 1597463007       ; 0x5f3759df
        sub     ecx, eax              ; 0x5f3759df - ...
        movd    xmm1, ecx
        mulss   xmm0, dword ptr [rip + .LCPI0_0]    ; x *= 0.5
        mulss   xmm0, xmm1                          ; x *= y
        mulss   xmm0, xmm1                          ; x *= y
        addss   xmm0, dword ptr [rip + .LCPI0_1]    ; x += 1.5
        mulss   xmm0, xmm1                          ; x *= y
        ret

我没有找到任何参考程序集（如果您有，请告诉我！），但它对我来说似乎相当不错。我只是不确定为什么浮点数被移动到eax中，然后进行移位和整数减法。也许SSE寄存器不支持这些操作？

使用-O3的clang 9.0将C代码编译为基本相同的汇编代码。所以这是一个好迹象。

---

值得指出的是，如果您真的想在实践中使用它：请不要。正如评论中所指出的，现代x86 CPU具有用于此功能的专用指令，比这个hack更快且更准确。不幸的是，1.0 / x.sqrt()似乎无法优化到该指令。因此，如果您真的需要速度，则使用_mm_rsqrt_ps内部函数可能是可行的方法。然而，这又需要使用unsafe代码。我不会在这个答案中详细介绍，因为只有少数程序员实际上需要它。

- Lukas Kalbertodt

---

根据英特尔指令集手册，没有只移动128位寄存器中最低的32位（类似于addss或mulss）的整数移位操作。但是如果可以忽略xmm0寄存器中其他96位，则可以使用psrld指令来实现。整数减法也是如此。 - fsasm

我承认对于 Rust 我几乎一无所知，但是 "unsafe" 不是快速反平方根的核心属性吗？它完全不尊重数据类型等。 - Gloweye

@Gloweye，我们讨论的“不安全”是不同类型的。一个快速的近似值远离最佳点，而另一个则随意处理未定义的行为。 - Deduplicator

@Gloweye：从数学上讲，“fast_inv_sqrt”函数的最后一部分只是一个牛顿迭代步骤，用于寻找更好的“inv_sqrt”近似值。这一部分并不会存在任何安全问题。诡计在于第一部分，它会找到一个良好的近似值。“fast_inv_sqrt”之所以有效，是因为它对浮点数的指数部分进行了整除2操作，事实上，“sqrt(pow(0.5,x))=pow(0.5,x/2)”。 - MSalters

@fsasm：没错；将movd到EAX再返回是当前编译器中的一个被忽略的优化。（而且，调用约定在XMM的低元素中传递/返回标量float并允许高位为垃圾。但请注意，如果它*被零扩展，它可以很容易地保持这种方式：右移不会引入非零元素，从_mm_set_epi32(0,0,0,0x5f3759df)减去也不会，即一个movd加载。在psrld之前需要一个movdqa xmm1,xmm0来复制寄存器。从FP指令转发到整数和反之的旁路延迟被mulss延迟隐藏了。 - Peter Cordes

Answer 2

这个例子是用 Rust 中不太常见的 union 实现的：

union FI {
    f: f32,
    i: i32,
}

fn inv_sqrt(x: f32) -> f32 {
    let mut u = FI { f: x };
    unsafe {
        u.i = 0x5f3759df - (u.i >> 1);
        u.f * (1.5 - 0.5 * x * u.f * u.f)
    }
}

---

在一个 x86-64 Linux 系统上，使用 criterion 库进行了微基准测试。令人惊讶的是，在 Rust 中自带的 sqrt().recip() 是最快的。但是，当然，任何微基准测试的结果都应该持怀疑态度。

---

inv sqrt with transmute time:   [1.6605 ns 1.6638 ns 1.6679 ns]
inv sqrt with union     time:   [1.6543 ns 1.6583 ns 1.6633 ns]
inv sqrt with to and from bits
                        time:   [1.7659 ns 1.7677 ns 1.7697 ns]
inv sqrt with powf      time:   [7.1037 ns 7.1125 ns 7.1223 ns]
inv sqrt with sqrt then recip
                        time:   [1.5466 ns 1.5488 ns 1.5513 ns]

- edwardw

---

我毫不惊讶 sqrt().inv() 是最快的。现在sqrt和inv都是单条指令，并且执行速度非常快。而Doom是在不能保证硬件浮点运算存在的时代编写的，像sqrt这样的超越函数肯定是通过软件实现的。赞成这个基准测试结果。 - Martin Bonner supports Monica

让我惊讶的是，transmute 显然与 to_ 和 from_bits 不同 -- 即使在优化之前，我也会认为它们是等效的指令。 - trent

@MartinBonner（另外，虽然这并不重要，但sqrt不是超越函数。） - benrg

任何支持除法的硬件浮点单元通常也支持平方根。IEEE“基本”操作（+ - * / sqrt）需要产生正确舍入的结果；这就是为什么SSE提供了所有这些操作，但不提供exp、sin或其他操作的原因。实际上，除法和平方根通常在同一个执行单元上运行，并且设计方式相似。请参见HW div/sqrt unit details。无论如何，它们与乘法相比仍然不快，特别是在延迟方面。 - Peter Cordes

无论如何，Skylake在除法/平方根方面的流水线比以前的uarches显着改善。请参见Floating point division vs floating point multiplication中Agner Fog表格的一些摘录。如果您在循环中没有做太多其他工作，因此sqrt + div是瓶颈，您可能希望使用HW快速倒数平方根（而不是quake hack）+牛顿迭代。特别是对于FMA，这对吞吐量非常有用，但对延迟则不然。Fast vectorized rsqrt and reciprocal with SSE/AVX depending on precision - Peter Cordes

Answer 3

您可以使用std::mem::transmute进行所需的转换：

fn inv_sqrt(x: f32) -> f32 {
    let xhalf = 0.5f32 * x;
    let mut i: i32 = unsafe { std::mem::transmute(x) };
    i = 0x5f3759df - (i >> 1);
    let mut res: f32 = unsafe { std::mem::transmute(i) };
    res = res * (1.5f32 - xhalf * res * res);
    res
}

您可以在此处查看实时示例：这里。

- Deedee Megadoodoo

---

没有安全问题是没有问题的，但有一种方法可以在不使用显式的unsafe块的情况下实现这一点，因此我建议使用 f32::to_bits 和 f32::from_bits 重写这个答案。与transmute不同，它清晰地表达了意图，而大多数人可能会认为transmute是“魔法”。 - user11877195

@Sahsahae，我刚刚发布了一个使用你提到的两个函数的答案 :) 我同意unsafe在这里是不必要的，应该避免使用。 - Lukas Kalbertodt