¿Es posible escribir la función rápida InvSqrt () de Quake en Rust?

Esto es solo para satisfacer mi propia curiosidad.

¿Hay una implementación de esto?

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;
}

en el óxido? Si existe, publique el código.

Lo intenté y fallé. No sé cómo codificar el número flotante con formato entero. Aquí está mi intento:

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
}

Referencia:
1. Origen de Quake3's Fast InvSqrt () - Página 1
2. Comprensión de Quake's Fast Inverse Square Root
3. FAST INVERSE SQUARE ROOT.pdf
4. código fuente: q_math.c # L552-L572

No sé cómo codificar el número flotante con formato entero.

Hay una función para eso: f32::to_bitsque devuelve un u32. También existe la función para la otra dirección: f32::from_bitsque toma un u32argumento. Estas funciones son preferibles a mem::transmutelas últimas, unsafey es difícil de usar.

Con eso, aquí está la implementación de 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)
}

( Área de juegos )

Esta función se compila en el siguiente ensamblado en 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

No he encontrado ningún ensamblaje de referencia (si lo ha hecho, ¡dígamelo!), Pero me parece bastante bueno. No estoy seguro de por qué se movió el flotador eaxsolo para hacer el cambio y la resta de enteros. ¿Quizás los registros SSE no admiten esas operaciones?

clang 9.0 con -O3compila el código C para básicamente el mismo ensamblado . Entonces esa es una buena señal.

Vale la pena señalar que si realmente quieres usar esto en la práctica: por favor no lo hagas. Como señaló benrg en los comentarios , las CPU modernas x86 tienen una instrucción especializada para esta función que es más rápida y precisa que este truco. Desafortunadamente, 1.0 / x.sqrt() no parece optimizar a esa instrucción . Entonces, si realmente necesita la velocidad, usar los _mm_rsqrt_psintrínsecos es probablemente el camino a seguir. Sin embargo, esto nuevamente requiere unsafecódigo. No entraré en muchos detalles en esta respuesta, ya que una minoría de programadores realmente lo necesitará.

Este se implementa con menos conocido unionen 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)
    }
}

Hice algunos micro benchmarks usando criterioncajón en una caja Linux x86-64. Sorprendentemente, el propio Rust sqrt().recip()es el más rápido. Pero, por supuesto, cualquier resultado micro referencial debe tomarse con un grano de sal.

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]

Puede usar std::mem::transmutepara hacer la conversión necesaria:

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
}

Puedes buscar un ejemplo en vivo aquí: aquí