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

101

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

Flyq
fuente
44
Según tengo entendido, este código es UB en C debido a la violación de la estricta regla de alias . La forma estándar de realizar este tipo de castigo es con a union.
trentcl
44
@trentcl: Yo tampoco creo que unionfuncione. memcpydefinitivamente funciona, aunque es detallado.
Matthieu M.
14
@MatthieuM. El tipo de punteo con uniones es perfectamente válido C , pero no válido C ++.
Moira
44
Supongo que esta pregunta está bien desde una perspectiva de pura curiosidad, pero comprenda que los tiempos han cambiado. En x86, las instrucciones rsqrtssy rsqrtps, introducidas con el Pentium III en 1999, son más rápidas y precisas que este código. BRAZO NEON tiene vrsqrteque es similar. Y cualesquiera que sean los cálculos que utilizó Quake III para esto, probablemente se realizarían en la GPU en estos días de todos modos.
benrg

Respuestas:

87

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á.

Lukas Kalbertodt
fuente
44
De acuerdo con la Guía de Intel intrínseco no hay ninguna operación de cambio de número entero que sólo desplaza el más bajo de 32 bits del registro analógico de 128 bits a addsso mulss. Pero si los otros 96 bits de xmm0 pueden ignorarse, entonces uno podría usar la psrldinstrucción. Lo mismo ocurre con la resta de enteros.
fsasm
¿Admitiré saber casi nada sobre el óxido, pero no es "inseguro" básicamente una propiedad central de fast_inv_sqrt? Con su total falta de respeto por los tipos de datos y demás.
Gloweye
12
@Gloweye Sin embargo, es un tipo diferente de "inseguro" del que hablamos. Una aproximación rápida que obtiene un mal valor demasiado lejos del punto óptimo, frente a algo que juega rápido y suelto con un comportamiento indefinido.
Deduplicador
8
@Gloweye: Matemáticamente, la última parte de eso fast_inv_sqrtes solo un paso de iteración de Newton-Raphson para encontrar una mejor aproximación de inv_sqrt. No hay nada inseguro en esa parte. El truco está en la primera parte, que encuentra una buena aproximación. Eso funciona porque está haciendo un número entero dividido por 2 en la parte exponente del flotador y, de hechosqrt(pow(0.5,x))=pow(0.5,x/2)
MSalters el
1
@fsasm: Eso es correcto; movdEAX y viceversa es una optimización perdida por los compiladores actuales. (Y sí, pasan convenciones de llamada / retorno escalar floaten el elemento de baja de un XMM y permiten altos bits a la basura Pero tenga en cuenta que si. Fue extendida a cero, puede mantenerse fácilmente de esa manera: desplazamiento de la derecha no introduce no cero elementos y tampoco resta de _mm_set_epi32(0,0,0,0x5f3759df), es decir, una movdcarga. Necesitaría un movdqa xmm1,xmm0para copiar el registro antes psrld. Evitar la latencia del reenvío de instrucciones FP a entero y viceversa está oculto por la mulsslatencia.
Peter Cordes
37

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]
edwardw
fuente
22
No estoy para nada sorprendido, sqrt().inv()es el más rápido. Tanto sqrt como inv son instrucciones únicas en estos días, y van bastante rápido. Doom se escribió en los días en que no era seguro asumir que había un punto flotante de hardware, y que las funciones trascendentales como sqrt definitivamente habrían sido software. +1 para los puntos de referencia.
Martin Bonner apoya a Mónica el
44
Lo que me sorprende es que transmuteaparentemente es diferente to_y from_bitsespero que sean equivalentes a las instrucciones incluso antes de la optimización.
trentcl
2
@MartinBonner (Además, no es que importe, pero sqrt no es una función trascendental .)
benrg
44
@ MartininBonner: cualquier FPU de hardware que admita división normalmente también admitirá sqrt. Las operaciones "básicas" IEEE (+ - * / sqrt) son necesarias para producir un resultado correctamente redondeado; es por eso que SSE proporciona todas esas operaciones pero no exp, sin, o lo que sea. De hecho, divide y sqrt generalmente se ejecutan en la misma unidad de ejecución, diseñada de manera similar. Vea los detalles de la unidad HW div / sqrt . De todos modos, todavía no son rápidos en comparación con la multiplicación, especialmente en latencia.
Peter Cordes el
1
De todos modos, Skylake tiene una canalización significativamente mejor para div / sqrt que los uarches anteriores. Ver División de punto flotante versus multiplicación de punto flotante para algunos extractos de la tabla de Agner Fog. Si no está haciendo mucho otro trabajo en un bucle, por lo que sqrt + div es un cuello de botella, es posible que desee usar HW recíproco rápido sqrt (en lugar del terremoto) + una iteración de Newton. Especialmente con FMA que es bueno para el rendimiento, si no latencia. Fast vectorizado rsqrt y recíproca con SSE / AVX dependiendo de la precisión
Peter Cordes
10

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í

Muy fresco
fuente
44
No hay nada de malo en inseguro, pero hay una manera de hacerlo sin un bloqueo explícito inseguro, por lo que sugeriría reescribir esta respuesta usando f32::to_bitsy f32::from_bits. También lleva la intención claramente diferente a transmutar, que la mayoría de la gente probablemente ve como "magia".
Sahsahae
55
@Sahsahae Acabo de publicar una respuesta usando las dos funciones que mencionaste :) Y estoy de acuerdo, unsafedebería evitarse aquí, ya que no es necesario.
Lukas Kalbertodt