是否可以在 Rust 中编写雷神之锤的快速 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;
}

生锈了吗？如果存在，发布代码。

我试过了，失败了。不知道如何用整数格式对浮点数进行编码。以下是我的尝试:

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.地震起源3的Fast InvSqrt（）-第1页
2.了解地震的快速平方反比根
3.快速反比平方根。pdf
4.源代码：q_math.c#L552-L572

Answer 1

您可以使用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
}

你可以在这里寻找一个活生生的例子:这里

Answer 2

这个是在 Rust 中用鲜为人知的并集实现的：

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盒子上使用标准板条箱做了一些微观基准测试。令人惊讶的是，鲁斯特自己的平方（）.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]

Answer 3

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

有一个函数：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寄存器不支持这些操作？

clang 9.0用< code>-O3把C代码编译成基本相同的汇编。所以这是个好迹象。

值得指出的是，如果你真的想在实践中使用这个：请不要。正如本格在评论中指出的那样，现代x86 CPU有一个专门的指令来实现这个功能，它比这个方法更快、更准确。遗憾的是，1.0/x.sqrt（）似乎没有针对该指令进行优化。因此，如果您真的需要速度，使用_mm_rsqrt_psintrinsic可能是一种方法。然而，这同样需要不安全代码。在这个答案中，我不会详细讨论，因为少数程序员实际上需要它。