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

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

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

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]

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