Ist es möglich, Quakes schnelle InvSqrt () -Funktion in Rust zu schreiben?

101

Dies ist nur, um meine eigene Neugier zu befriedigen.

Gibt es eine Implementierung davon:

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

in Rust? Wenn es existiert, geben Sie den Code ein.

Ich habe es versucht und bin gescheitert. Ich weiß nicht, wie ich die Gleitkommazahl im Ganzzahlformat codieren soll. Hier ist mein Versuch:

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
}

Referenz:
1. Ursprung von Quake3s Fast InvSqrt () - Page 1
2. Grundlegendes zu Quakes Fast Inverse Square Root
3. FAST INVERSE SQUARE ROOT.pdf
4. Quellcode: q_math.c # L552-L572

Flyq
quelle
die C # -Version
Flyq
4
Soweit ich weiß, ist dieser Code UB in C, da er gegen die strikte Aliasing-Regel verstößt . Die Standard-gesegnete Art, diese Art von Punning durchzuführen, ist mit a union.
Trentcl
4
@trentcl: Ich denke auch nicht unionfunktioniert. memcpyfunktioniert definitiv, obwohl es ausführlich ist.
Matthieu M.
14
@MatthieuM. Typ Punning mit Gewerkschaften ist vollkommen gültiges C , aber nicht gültiges C ++.
Moira
4
Ich nehme an, diese Frage ist aus rein neugieriger Sicht in Ordnung, aber bitte haben Sie Verständnis dafür, dass sich die Zeiten geändert haben. Auf x86, die rsqrtssund rsqrtpsAnweisungen, im Jahr 1999 mit dem Pentium III eingeführt, sind schneller und genauer als dieser Code. ARM NEON hat vrsqrtewas ähnlich ist. Und für welche Berechnungen Quake III dies auch verwendete, würde es heutzutage wahrscheinlich sowieso auf der GPU durchgeführt werden.
Benrg

Antworten:

87

Ich weiß nicht, wie ich die Gleitkommazahl im Ganzzahlformat codieren soll.

Dafür gibt es eine Funktion: f32::to_bitsdie ein zurückgibt u32. Es gibt auch die Funktion für die andere Richtung: f32::from_bitsdie ein u32als Argument nimmt. Diese Funktionen werden gegenüber mem::transmuteletzteren bevorzugt unsafeund sind schwierig zu bedienen.

Damit ist hier die Implementierung von 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)
}

( Spielplatz )


Diese Funktion wird auf x86-64 in die folgende Assembly kompiliert:

.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

Ich habe keine Referenzbaugruppe gefunden (wenn ja, bitte sagen Sie es mir!), Aber es scheint mir ziemlich gut zu sein. Ich bin mir nur nicht sicher, warum der Float nur bewegt wurde eax, um die Verschiebung und die Ganzzahlsubtraktion durchzuführen. Vielleicht unterstützen SSE-Register diese Operationen nicht?

clang 9.0 mit -O3kompiliert den C-Code auf im Grunde dieselbe Assembly . Das ist also ein gutes Zeichen.


Es sei darauf hingewiesen, dass Sie dies bitte nicht tun sollten, wenn Sie dies tatsächlich in der Praxis anwenden möchten. Wie benrg in den Kommentaren hervorhob , verfügen moderne x86-CPUs über eine spezielle Anweisung für diese Funktion, die schneller und genauer als dieser Hack ist. Leider 1.0 / x.sqrt() scheint sich diese Anweisung nicht zu optimieren . Wenn Sie also wirklich die Geschwindigkeit brauchen, ist die Verwendung der _mm_rsqrt_psIntrinsics wahrscheinlich der richtige Weg. Dies erfordert jedoch wieder unsafeCode. Ich werde in dieser Antwort nicht ins Detail gehen, da eine Minderheit der Programmierer sie tatsächlich brauchen wird.

Lukas Kalbertodt
quelle
4
Laut Intel Intrinsics Guide gibt es keine Ganzzahlverschiebungsoperation, die nur das niedrigste 32-Bit des 128-Bit-Registers analog zu addssoder verschiebt mulss. Aber wenn die anderen 96 Bits von xmm0 ignoriert werden können, könnte man den psrldBefehl verwenden. Gleiches gilt für die Ganzzahlsubtraktion.
Fsasm
Ich gebe zu, so gut wie nichts über Rost zu wissen, aber ist "unsicher" nicht im Grunde eine Kerneigenschaft von fast_inv_sqrt? Mit seiner völligen Missachtung von Datentypen und dergleichen.
Gloweye
12
@Gloweye Es ist eine andere Art von "unsicher", über die wir jedoch sprechen. Eine schnelle Annäherung, die einen schlechten Wert zu weit vom Sweet Spot entfernt hat, im Gegensatz zu etwas, das schnell und locker mit undefiniertem Verhalten spielt.
Deduplikator
8
@Gloweye: Mathematisch gesehen ist der letzte Teil davon fast_inv_sqrtnur ein Newton-Raphson-Iterationsschritt, um eine bessere Annäherung an zu finden inv_sqrt. An diesem Teil ist nichts Unsicheres. Der Trick ist im ersten Teil, der eine gute Annäherung findet. Das funktioniert, weil es eine ganzzahlige Division durch 2 auf dem Exponententeil des sqrt(pow(0.5,x))=pow(0.5,x/2)
Floats macht
1
@fsasm: Das stimmt; movdzu EAX und zurück ist eine verpasste Optimierung durch aktuelle Compiler. (Und ja, das Aufrufen von Konventionen übergibt / gibt den Skalar floatim Low-Element eines XMM durch und lässt zu, dass High-Bits Müll sind. Beachten Sie jedoch , dass dies leicht so bleiben kann , wenn es auf Null erweitert wurde: Rechtsverschiebung führt nicht zu Nicht-Verschiebungen. Null Elemente und auch keine Subtraktion von _mm_set_epi32(0,0,0,0x5f3759df), dh eine movdLast. Sie müssten movdqa xmm1,xmm0die Registrierung vorher kopieren psrld. Die Umgehungslatenz von der Weiterleitung des FP-Befehls zur Ganzzahl und umgekehrt wird durch die mulssLatenz verborgen .
Peter Cordes
37

Dieser ist mit weniger bekannt unionin Rust implementiert :

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

Habe einige Micro-Benchmarks mit criterioncrate auf einer x86-64-Linux-Box durchgeführt. Überraschenderweise ist Rusts eigene sqrt().recip()die schnellste. Aber natürlich sollte jedes Mikro-Benchmark-Ergebnis mit einem Körnchen Salz aufgenommen werden.

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
quelle
22
Ich bin nicht im geringsten überrascht, sqrt().inv()ist am schnellsten. Sowohl sqrt als auch inv sind heutzutage einzelne Anweisungen und gehen ziemlich schnell. Doom wurde in den Tagen geschrieben, als es nicht sicher war anzunehmen, dass es überhaupt Hardware-Gleitkomma gab, und transzendentale Funktionen wie sqrt wären definitiv Software gewesen. +1 für die Benchmarks.
Martin Bonner unterstützt Monica
4
Was mich überrascht, ist, dass transmutees sich anscheinend von to_und unterscheidet from_bits- ich würde erwarten, dass diese bereits vor der Optimierung den Anweisungen entsprechen.
Trentcl
2
@ MartinBonner (Auch nicht, dass es wichtig ist, aber sqrt ist keine transzendentale Funktion .)
Benrg
4
@MartinBonner: Jede Hardware-FPU, die die Division unterstützt, unterstützt normalerweise auch sqrt. IEEE "grundlegende" Operationen (+ - * / sqrt) sind erforderlich, um ein korrekt gerundetes Ergebnis zu erzielen. Deshalb bietet SSE all diese Operationen an, aber nicht exp, sin oder was auch immer. Tatsächlich werden dividieren und sqrt normalerweise auf derselben Ausführungseinheit ausgeführt, die auf ähnliche Weise entworfen wurde. Siehe Details zur HW div / sqrt-Einheit . Auf jeden Fall sind sie im Vergleich zur Multiplikation immer noch nicht schnell, insbesondere in Bezug auf die Latenz.
Peter Cordes
1
Wie auch immer, Skylake hat ein deutlich besseres Pipelining für div / sqrt als frühere Uarches. Einige Auszüge aus der Tabelle von Agner Fog finden Sie unter Gleitkommadivision vs. Gleitkommamultiplikation . Wenn Sie nicht viel andere Arbeit in einer Schleife erledigen, sodass sqrt + div ein Engpass ist, möchten Sie möglicherweise HW schnelles reziprokes sqrt (anstelle des Quake-Hack) + eine Newton-Iteration verwenden. Besonders bei FMA ist das gut für den Durchsatz, wenn nicht für die Latenz. Schnelle vektorisierte rsqrt und wechselseitig mit SSE / AVX je nach Präzision
Peter Cordes
10

Sie können verwenden std::mem::transmute, um die erforderliche Konvertierung vorzunehmen:

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
}

Hier können Sie nach einem Live-Beispiel suchen: hier

Echt frisch
quelle
4
Es ist nichts Falsches an unsicher, aber es gibt eine Möglichkeit, dies ohne explizite unsichere Blockierung zu tun. Daher würde ich vorschlagen, diese Antwort mit f32::to_bitsund neu zu schreiben f32::from_bits. Es trägt auch die Absicht eindeutig im Gegensatz zur Umwandlung, die die meisten Menschen wahrscheinlich als "Magie" betrachten.
Sahsahae
5
@Sahsahae Ich habe gerade eine Antwort mit den beiden von Ihnen erwähnten Funktionen gepostet :) Und ich stimme zu, unsafesollte hier vermieden werden, da es nicht notwendig ist.
Lukas Kalbertodt