rust - 是否可以用 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 中？如果存在，请发布代码。

我尝试了一下，失败了。我不知道如何使用整数格式对 float 进行编码。这是我的尝试:

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.Origin of Quake3's Fast InvSqrt() - Page 1
2.Understanding Quake’s Fast Inverse Square Root
3.FAST INVERSE SQUARE ROOT.pdf
4.source code: q_math.c#L552-L572

Answer 1

I don't know how to encode the float number using integer format.

有一个函数可以实现这一点: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)
}

(Playground)

<小时/>

此函数在 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 与 -O3 将 C 代码编译为 basically the same assembly 。所以这是一个好兆头。

<小时/>

值得指出的是，如果您确实想在实践中使用它:请不要这样做。作为 benrg pointed out in the comments ，现代 x86 CPU 有一个专门用于此功能的指令，它比此 hack 更快、更准确。不幸的是，1.0/x.sqrt() does not seem to optimize to that instruction 。因此，如果您确实需要速度，那么使用 the _mm_rsqrt_ps intrinsics 可能是最佳选择。然而，这又需要不安全代码。我不会在这个答案中详细介绍，因为少数程序员实际上需要它。