Est-il possible d'écrire la fonction rapide InvSqrt () de Quake dans Rust?

101

C'est juste pour satisfaire ma propre curiosité.

Y a-t-il une implémentation de ceci:

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? S'il existe, affichez le code.

J'ai essayé et j'ai échoué. Je ne sais pas comment encoder le nombre flottant en utilisant le format entier. Voici ma tentative:

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
}

Référence:
1. Origine de Quake3's Fast InvSqrt () - Page 1
2. Comprendre la racine carrée inverse rapide de Quake
3. FAST INVERSE SQUARE ROOT.pdf
4. code source: q_math.c # L552-L572

Flyq
la source
4
Si je comprends bien, ce code est UB en C en raison de la violation de la règle stricte d'alias . La manière standard bénie d'effectuer ce type de punition de type est avec a union.
trentcl
4
@trentcl: Je ne pense pas que ça unionmarche non plus. memcpyfonctionne certainement, bien qu'il soit verbeux.
Matthieu M.
14
@MatthieuM. Le type punning avec unions est parfaitement C valide , mais pas valide C ++.
Moira
4
Je suppose que cette question est correcte du point de vue de la pure curiosité, mais veuillez comprendre que les temps ont changé. Sur x86, les instructions rsqrtsset rsqrtps, introduites avec le Pentium III en 1999, sont plus rapides et plus précises que ce code. ARM NEON a vrsqrtequi est similaire. Et quels que soient les calculs utilisés par Quake III, cela serait probablement fait sur le GPU ces jours-ci de toute façon.
benrg

Réponses:

87

Je ne sais pas comment encoder le nombre flottant en utilisant le format entier.

Il y a une fonction pour cela: f32::to_bitsqui retourne un u32. Il y a aussi la fonction pour l'autre direction: f32::from_bitsqui prend un u32argument comme. Ces fonctions sont préférées par rapport mem::transmuteà cette dernière unsafeet délicate à utiliser.

Avec cela, voici la mise en œuvre 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)
}

( Aire de jeux )


Cette fonction se compile vers l'assembly suivant sur 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

Je n'ai trouvé aucun assemblage de référence (si vous en avez, dites-le moi!), Mais il me semble assez bon. Je ne sais tout simplement pas pourquoi le flotteur a été déplacé eaxjuste pour effectuer le décalage et la soustraction d'entiers. Peut-être que les registres SSE ne prennent pas en charge ces opérations?

clang 9.0 avec -O3compile le code C en gros dans le même assemblage . C'est donc un bon signe.


Il convient de souligner que si vous souhaitez réellement l'utiliser dans la pratique: veuillez ne pas le faire. Comme l'a souligné benrg dans les commentaires , les processeurs x86 modernes ont une instruction spécialisée pour cette fonction qui est plus rapide et plus précise que ce hack. Malheureusement, 1.0 / x.sqrt() ne semble pas optimiser cette instruction . Donc, si vous avez vraiment besoin de vitesse, l'utilisation de l' _mm_rsqrt_psintrinsèque est probablement la voie à suivre. Cependant, cela nécessite à nouveau du unsafecode. Je n'entrerai pas dans les détails de cette réponse, car une minorité de programmeurs en aura réellement besoin.

Lukas Kalbertodt
la source
4
Selon le Guide Intel Intrinsics, il n'y a pas d'opération de décalage d'entier qui décale uniquement les 32 bits les plus bas du registre analogique de 128 bits vers addssou mulss. Mais si les 96 autres bits de xmm0 peuvent être ignorés, alors on pourrait utiliser l' psrldinstruction. Il en va de même pour la soustraction entière.
fsasm
J'avoue ne rien savoir de la rouille, mais n'est-ce pas "dangereux" fondamentalement une propriété de base de fast_inv_sqrt? Avec son manque total de respect pour les types de données et autres.
Gloweye
12
@Gloweye C'est un type différent de "dangereux" dont nous parlons cependant. Une approximation rapide qui obtient une mauvaise valeur trop loin du sweet spot, par rapport à un jeu rapide et lâche avec un comportement indéfini.
Déduplicateur
8
@Gloweye: Mathématiquement, la dernière partie de cela fast_inv_sqrtn'est qu'une étape d'itération de Newton-Raphson pour trouver une meilleure approximation de inv_sqrt. Il n'y a rien de dangereux dans cette partie. La ruse est dans la première partie, qui trouve une bonne approximation. Cela fonctionne parce qu'il fait une division entière par 2 sur la partie exposante du flotteur, et en effetsqrt(pow(0.5,x))=pow(0.5,x/2)
MSalters
1
@fsasm: C'est correct; movdà EAX et retour est une optimisation manquée par les compilateurs actuels. (Et oui, les conventions d'appel passent / retournent un scalaire floatdans l'élément bas d'un XMM et permettent aux bits élevés d'être des ordures. Mais notez que s'il était étendu à zéro, il peut facilement rester ainsi: le décalage à droite n'introduit zéro éléments et ni ne soustraction _mm_set_epi32(0,0,0,0x5f3759df), soit une movdcharge Vous auriez besoin d' un. movdqa xmm1,xmm0copier le reg avant psrldBypass temps d' attente de transfert d'instruction FP à l' entier et vice versa est caché par. mulsstemps d' attente.
Peter Cordes
37

Celui-ci est implémenté avec moins connu uniondans 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)
    }
}

A fait quelques micro benchmarks en utilisant une criterioncaisse sur une boîte Linux x86-64. Étonnamment Rust'ssqrt().recip() est le plus rapide. Mais bien sûr, tout résultat de micro-benchmark doit être pris avec un grain de sel.

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
la source
22
Je ne suis pas du tout surpris, sqrt().inv()c'est le plus rapide. Sqrt et inv sont des instructions uniques de nos jours, et vont assez vite. Doom a été écrit à l'époque où il n'était pas sûr de supposer qu'il y avait du matériel flottant du tout, et les fonctions transcendantales comme sqrt auraient certainement été des logiciels. +1 pour les benchmarks.
Martin Bonner soutient Monica
4
Ce qui me surprend, c'est qu'il transmuteest apparemment différent de to_et from_bits- je m'attends à ce que ceux-ci soient équivalents à l'instruction avant même l'optimisation.
trentcl
2
@MartinBonner (De plus, ce n'est pas important, mais sqrt n'est pas un fonction transcendantale .)
benrg
4
@MartinBonner: Tout FPU matériel qui prend en charge la division prend normalement également en charge sqrt. Les opérations "de base" IEEE (+ - * / sqrt) sont nécessaires pour produire un résultat correctement arrondi; c'est pourquoi SSE fournit toutes ces opérations mais pas exp, sin ou quoi que ce soit. En fait, divide et sqrt s'exécutent généralement sur la même unité d'exécution, conçue de la même manière. Voir les détails de l'unité HW div / sqrt . Quoi qu'il en soit, ils ne sont toujours pas rapides par rapport à la multiplication, en particulier en latence.
Peter Cordes
1
Quoi qu'il en soit, Skylake a un pipeline bien meilleur pour div / sqrt que les uarches précédents. Voir Division en virgule flottante vs multiplication en virgule flottante pour certains extraits du tableau d'Agner Fog. Si vous ne faites pas beaucoup d'autres travaux dans une boucle, donc sqrt + div est un goulot d'étranglement, vous voudrez peut-être utiliser sqrt réciproque rapide HW (au lieu du hack de tremblement de terre) + une itération de Newton. Surtout avec FMA qui est bon pour le débit, sinon la latence. Rsqrt vectorisé rapide et réciproque avec SSE / AVX selon la précision
Peter Cordes
10

Vous pouvez utiliser std::mem::transmutepour effectuer la conversion nécessaire:

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
}

Vous pouvez rechercher un exemple en direct ici: ici

Vraiment frais
la source
4
Il n'y a rien de mal avec dangereux, mais il y a un moyen de le faire sans bloc dangereux explicite, donc je suggère de réécrire cette réponse en utilisant f32::to_bitset f32::from_bits. Il porte également clairement l'intention contrairement à la transmutation, que la plupart des gens considèrent probablement comme «magique».
Sahsahae
5
@Sahsahae Je viens de poster une réponse en utilisant les deux fonctions que vous avez mentionnées :) Et je suis d'accord, unsafedevrait être évité ici, car ce n'est pas nécessaire.
Lukas Kalbertodt