Czy w Rust można napisać szybką funkcję InvSqrt () Quake'a?

101

Ma to zaspokoić moją ciekawość.

Czy istnieje implementacja tego:

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

w Rust? Jeśli istnieje, opublikuj kod.

Próbowałem i nie udało mi się. Nie wiem, jak zakodować liczbę zmiennoprzecinkową przy użyciu formatu liczb całkowitych. Oto moja próba:

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
}

Odniesienia:
1. Pochodzenie Quake3's Fast InvSqrt () - Strona 1
2. Zrozumienie Quake's Fast Inverse Square Root
3. FAST INVERSE SQUARE ROOT.pdf
4. Kod źródłowy: q_math.c # L552-L572

Flyq
źródło
4
Jak rozumiem, ten kod jest UB w C z powodu naruszenia ścisłej zasady aliasingu . Standardowo pobłogosławionym sposobem wykonywania tego rodzaju znakowania jest użycie znaku union.
trentcl
4
@trentcl: Też nie sądzę, że uniondziała. memcpyzdecydowanie działa, choć jest to pełne.
Matthieu M.
14
@MatthieuM. Pisanie na znakach za pomocą związków jest całkowicie poprawnym C , ale nie poprawnym C ++.
Moira
4
Przypuszczam, że to pytanie jest w porządku z czystej ciekawości, ale proszę zrozumieć, że czasy się zmieniły. Na x86, to rsqrtssi rsqrtpsinstrukcje, wprowadzone z Pentium III w 1999 roku, są szybsze i dokładniejsze niż tego kodu. ARM NEON ma vrsqrteto, co jest podobne. I niezależnie od tego, jakie obliczenia zastosował Quake III, prawdopodobnie i tak zostałyby wykonane na GPU.
benrg

Odpowiedzi:

87

Nie wiem, jak zakodować liczbę zmiennoprzecinkową przy użyciu formatu liczb całkowitych.

Jest na to funkcja: f32::to_bitsktóra zwraca an u32. Istnieje również funkcja dla drugiego kierunku: f32::from_bitsktóry przyjmuje u32argument jako argument. Te funkcje są lepsze niż mem::transmutete ostatnie unsafei są trudne w użyciu.

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

( Plac zabaw )


Ta funkcja kompiluje się do następującego zestawu na 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

Nie znalazłem żadnego zestawu referencyjnego (jeśli tak, proszę powiedz mi!), Ale wydaje mi się, że jest całkiem dobry. Nie jestem tylko pewien, dlaczego zmiennoprzecinkowe zostało przeniesione eaxtylko po to, aby wykonać przesunięcie i odejmowanie liczb całkowitych. Może rejestry SSE nie obsługują tych operacji?

clang 9.0 z -O3kompiluje kod C w zasadzie do tego samego zestawu . To dobry znak.


Warto zauważyć, że jeśli rzeczywiście chcesz to wykorzystać w praktyce: nie rób tego. Jak zauważył Benrg w komentarzach , nowoczesne procesory x86 mają specjalną instrukcję dla tej funkcji, która jest szybsza i dokładniejsza niż ten hack. Niestety 1.0 / x.sqrt() nie wydaje się, aby optymalizować tę instrukcję . Jeśli więc naprawdę potrzebujesz prędkości, prawdopodobnie skorzystaj z funkcji _mm_rsqrt_pswewnętrznych . To jednak wymaga jeszcze unsafekodu. Nie będę szczegółowo omawiał tej odpowiedzi, ponieważ mniejszość programistów faktycznie będzie jej potrzebować.

Lukas Kalbertodt
źródło
4
Według Intel Intrinsics Guide nie ma operacji przesunięcia liczb całkowitych, która przesuwa tylko najniższy 32-bit ze 128-bitowego rejestru analogowego na addsslub mulss. Ale jeśli pozostałe 96 bitów xmm0 można zignorować, można użyć psrldinstrukcji. To samo dotyczy odejmowania liczb całkowitych.
fsasm
Przyznaję, że prawie nic nie wiem o rdzy, ale czy „niebezpieczne” nie jest w zasadzie podstawową właściwością fast_inv_sqrt? Z całkowitym brakiem szacunku dla typów danych i tym podobnych.
Gloweye,
12
@Gloweye To inny rodzaj „niebezpiecznych”, o których mówimy. Szybkie przybliżenie, które odbiera złą wartość zbyt daleko od najsłabszego miejsca, w porównaniu do czegoś, co gra szybko i luźno z nieokreślonym zachowaniem.
Deduplicator
8
@ Gloweye: Matematycznie ostatnia część fast_inv_sqrtto tylko jeden krok iteracji Newtona-Raphsona, aby znaleźć lepsze przybliżenie inv_sqrt. W tej części nie ma nic niebezpiecznego. Sztuczka znajduje się w pierwszej części, która znajduje dobre przybliżenie. To działa, ponieważ wykonuje dzielenie przez liczbę całkowitą przez 2 w części wykładniczej liczby zmiennoprzecinkowej, i rzeczywiściesqrt(pow(0.5,x))=pow(0.5,x/2)
MSalters
1
@fsasm: Zgadza się; movddo EAX iz powrotem jest brakującą optymalizacją obecnych kompilatorów. (I tak, wywoływanie konwencji przekazuje / zwraca skalar floatw dolnym elemencie XMM i pozwala na wyrzucanie dużych bitów. Pamiętaj jednak, że jeśli był rozszerzony do zera, może z łatwością pozostać w ten sposób: prawe przesunięcie nie wprowadza zero elementów i żadne nie odejmuje _mm_set_epi32(0,0,0,0x5f3759df), tj movd. obciążenia. Trzeba by wcześniej movdqa xmm1,xmm0skopiować reg psrld. Pominięcie opóźnienia z przekazania instrukcji FP do liczby całkowitej i odwrotnie jest ukryte przez mulssopóźnienie
Peter Cordes
37

Ten jest zaimplementowany z mniej znanym unionw 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)
    }
}

Czy niektóre mikro testy porównawcze przy użyciu criterionskrzynki na komputerze z systemem Linux x86-64. Zaskakująco własny sqrt().recip()jest najszybszy. Ale oczywiście każdy wynik mikroprocesora powinien być wzięty z odrobiną soli.

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
źródło
22
Nie jestem zaskoczony, że sqrt().inv()jest najszybszy. Zarówno sqrt, jak i inv są obecnie pojedynczymi instrukcjami i działają dość szybko. Doom został napisany w czasach, gdy nie było bezpiecznie zakładać, że w ogóle występuje zmiennoprzecinkowy sprzęt, a funkcje transcendentalne, takie jak sqrt, zdecydowanie byłyby oprogramowaniem. +1 za testy porównawcze.
Martin Bonner wspiera Monikę
4
Co zaskakuje mnie to, że transmutejest zupełnie różne od działania to_i from_bits- Spodziewam się tych instrukcji równoważne nawet przed optymalizacji.
trentcl
2
@MartinBonner (Poza tym nie ma to znaczenia, ale sqrt nie jest funkcją transcendentalną ).
benrg
4
@MartinBonner: Każdy sprzętowy układ FPU obsługujący podział zwykle obsługuje również sqrt. „Podstawowe” operacje IEEE (+ - * / sqrt) są wymagane do uzyskania poprawnie zaokrąglonego wyniku; dlatego SSE zapewnia wszystkie te operacje, ale nie exp, grzech lub cokolwiek innego. W rzeczywistości divide i sqrt zwykle działają na tej samej jednostce wykonawczej, zaprojektowanej w podobny sposób. Zobacz szczegółowe informacje na temat jednostki div / sqrt HW . W każdym razie nadal nie są szybkie w porównaniu do mnożenia, szczególnie w przypadku opóźnienia.
Peter Cordes,
1
W każdym razie Skylake ma znacznie lepszy potok dla div / sqrt niż poprzednie uarche. Zobacz dzielenie zmiennoprzecinkowe vs mnożenie zmiennoprzecinkowe, aby zapoznać się z niektórymi fragmentami tabeli Agner Fog. Jeśli nie wykonujesz zbyt wiele innych prac w pętli, więc sqrt + div jest wąskim gardłem, możesz użyć HW szybkiego odwrotnego sqrt (zamiast hakowania) + iteracji Newtona. Zwłaszcza z FMA, która jest dobra pod względem przepustowości, jeśli nie opóźnienia. Szybki wektoryzowany rsqrt i wzajemność z SSE / AVX w zależności od precyzji
Peter Cordes
10

Możesz użyć std::mem::transmutedo dokonania niezbędnej konwersji:

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
}

Przykład na żywo możesz znaleźć tutaj: tutaj

Prawdziwie świeże
źródło
4
Nie ma nic złego w niebezpiecznym, ale jest sposób, aby to zrobić bez wyraźnego niebezpiecznego bloku, więc sugeruję przepisanie tej odpowiedzi za pomocą f32::to_bitsi f32::from_bits. Niesie też intencje wyraźnie odmienne od transmutacji, które większość ludzi prawdopodobnie uważa za „magię”.
Sahsahae,
5
@ Sahahahae Właśnie opublikowałem odpowiedź za pomocą dwóch wymienionych przez ciebie funkcji :) I zgadzam się, że unsafenależy tego unikać, ponieważ nie jest to konieczne.
Lukas Kalbertodt