2009-06-16
1.f/sqrtf(float): FPU vs Q3 hack vs SSE
Na wczorajszym seminarium (piątkowym w poniedziałek) pojawiło się pytanie: o ile szybszy jest szybki odwrotny pierwiastek kwadratowy z Quake 3 od normalnego 1.f/sqrtf(float)? Generalnie hack jest bardzo szczegółowo opisany, przedyskutowany i świetnie zbadany - warto przejrzeć papiery w przypisach do artykułu na wikipedii (szczególnie 4 i 5).
Ale nie wiedziałem jakie są konkretne liczby, dlatego napisałem prosty test mający to zmierzyć. Oprócz porównania ze zwykłym liczeniem na FPU postanowiłem to jeszcze porównać z taką samą szybką metodą, ale z dodatkowym przybliżeniem oraz z liczeniem na SSE (instrukcja rsqrtps - Packed Single-Precision Floating-Point Square Root Reciprocal). Wszystkie te metody są bardzo przybliżone - na FPU jest najdokładniej; mimo iż liczby są na floatach to na FPU wszystko jest liczone na rejestrach 80-bitowych. Za to na SSE wszystko jest liczone na liczbach 32-bitowych - gdzieś czytałem, że jest to realizowane podobną metodą co właśnie szybka odwrotny pierwiastek kwadratowy. Konkretne rezultaty poniżej.
Kilka uwag:
Jeżeli chodzi o metodę z Quake 3 to biorę minimalnie poprawioną wersję z magiczną stałą 0x5f375a86 wykorzystując wynik z Fast Inverse Square Root (pdf) - w tekście dużo detali oraz odpowiednia stała dla wersji na double. (W Understanding and improving fast inverse square root hack uzyskano prawie taką samą stałą: 0x5f375a84).
Niżej metodę na FPU nazwałem "Normal", a metodę szybką "Fast". Oprócz SSE jest jeszcze metoda szybka z dodatkową iteracją przybliżenia pod koniec, która nazywam "FastExtra" - na podstawie: Magical Square Root Implementation In Quake III .
Uwaga: flaga -O2 w gcc psuje kod szybkiej metody, dlatego poniższe testy kompilowałem z -O1.
TEST: 1.f/sqrtf(float)
Test wygląda tak: Biorę 1048576 pseudolosowych liczb z zakresu (0.0, 1000.0) (zawsze te same), które zapisuję w tablicy i 200 razy wykonuję obliczenie odwrotności pierwiastka kwadratowego na każdej z kolejnych liczb - czyli wykonuję 209715200 takich operacji. Wynik każdej operacji zapisuję w innej tablicy niż tablica wejściowa. Testuję kolejno każdą z metod. Czas jest liczony w sekundach. Błąd jest liczony w porównaniu z metodą na FPU.
Wyniki na laptopie z Pentium 4-M 2GHz, DDR-266:
Wyniki na desktopie z AMD Athlon X2 BE-2300, DDR2-800:
Wersja SSE jest wyraźnie najszybsza, ale największy błąd wyniósł nawet 3e-3. Warto zwrócić uwagę na to jak różny jest średni błąd na różnych maszynach, ale to chyba normalne przy SSE - będzie to wyraźniejsze w drugim teście.
TEST: Normalize(float3)
Ponieważ odwrotność pierwiastka wykorzystuję przede wszystkim do normalizowania trójwymiarowych wektorów, postanowiłem jeszcze przetestować jak każda z tych metod się do tego nadaje. Tutaj obliczając błąd sprawdzam jak daleka jest długość wyjściowych wektorów od pożądanej jedynki.
Wyniki na Pentium 4-M 2GHz, DDR-266:
Wyniki na AMD Athlon X2 BE-2300, z DDR2-800:
Przyspieszenie w przypadku SSE jest zadowalające - na szybszej maszynie wręcz bardzo. Na laptopie (Pentium 4M) wersja Fast jest całkiem blisko SSE, chociaż znowu ma wyraźnie gorszą dokładność. Jednak tam zyski z SSE najprawdopodobniej giną przesłonięte przez opóźnienia wywołane fatalnie wolną pamięcią.
Jak to zwykle przy SSE, ważne jest odpowiednie ułożenie danych tak by dało się obliczać cztery wektory na raz, co w niektórych warunkach może być kłopotliwe, jeżeli nie niemożliwe. Dla wszystkich metod ułożyłem dane wektor po wektorze - tzn. xyzxyzxy..., za to dla SSE dane są wcześniej przekształcone w ciąg: xxxxyyyyzzzzxxx... . We wszystkich przypadkach pamięć jest wyrównana do granicy 16-bajtowej. Jeżeli idzie o implementację dla SSE, to wygląda następująco:
Mięso w pętli jest tłumaczone przez gcc na:
Na koniec: 1.f/sqrtf(0.f)
Przy okazji normalizacji interesującym było jeszcze obejrzeć jak się te metody zachowują dla wektorów zerowych. FPU i SSE jako wynik 1/sqrt(0) zwracają inf (FPU bardzo powoli - można to jakoś przełączyć?), co dalej przemnożone przez 0.f daje nan. Za to w przypadku metod Fast i FastExtra dostałem bardzo duże wartości (Fast: 1.98e+19, FastExtra: 2.97e+19), które śmiało można przemnożyć przez 0. Wynik będzie taki, że normalizacja wektora na FPU i SSE daje wektor [NaN, NaN, NaN], za to pozostałe przybliżone metody dadzą [0.f, 0.f, 0.f]. Zależy od zastosowania, który wynik jest bardziej pożądany, chociaż pewno w większości przypadków lepiej jest standardowo dorzucić ifa po obliczeniu sumy kwadratów współrzędnych, ale może to skomplikować obliczenia na SSE.
Błądy ujemny i dodatnie:
Sprawdziłem jeszcze (na Pentium 4M), dla normalizacji wektora, jak dokładnie wyglądają najgorsze błędy dodatnie i ujemne:
Jak widać błąd ujemny jest minimalny - nawet mniejszy niż przy FPU. Błąd ujemny wskazuje na to o ile długość znormalizowanego wektora przekracza jedynkę. Czyli w przypadku metod Fast i FastExtra większość błędów powoduje wyprodukowanie wektorów trochę krótszych niż wektory normalne. Przy tak dużym błędzie absolutnym jest to fantastyczna własność - przynajmniej w moich zastosowaniach.
Ale nie wiedziałem jakie są konkretne liczby, dlatego napisałem prosty test mający to zmierzyć. Oprócz porównania ze zwykłym liczeniem na FPU postanowiłem to jeszcze porównać z taką samą szybką metodą, ale z dodatkowym przybliżeniem oraz z liczeniem na SSE (instrukcja rsqrtps - Packed Single-Precision Floating-Point Square Root Reciprocal). Wszystkie te metody są bardzo przybliżone - na FPU jest najdokładniej; mimo iż liczby są na floatach to na FPU wszystko jest liczone na rejestrach 80-bitowych. Za to na SSE wszystko jest liczone na liczbach 32-bitowych - gdzieś czytałem, że jest to realizowane podobną metodą co właśnie szybka odwrotny pierwiastek kwadratowy. Konkretne rezultaty poniżej.
Kilka uwag:
Jeżeli chodzi o metodę z Quake 3 to biorę minimalnie poprawioną wersję z magiczną stałą 0x5f375a86 wykorzystując wynik z Fast Inverse Square Root (pdf) - w tekście dużo detali oraz odpowiednia stała dla wersji na double. (W Understanding and improving fast inverse square root hack uzyskano prawie taką samą stałą: 0x5f375a84).
Niżej metodę na FPU nazwałem "Normal", a metodę szybką "Fast". Oprócz SSE jest jeszcze metoda szybka z dodatkową iteracją przybliżenia pod koniec, która nazywam "FastExtra" - na podstawie: Magical Square Root Implementation In Quake III .
Uwaga: flaga -O2 w gcc psuje kod szybkiej metody, dlatego poniższe testy kompilowałem z -O1.
TEST: 1.f/sqrtf(float)
Test wygląda tak: Biorę 1048576 pseudolosowych liczb z zakresu (0.0, 1000.0) (zawsze te same), które zapisuję w tablicy i 200 razy wykonuję obliczenie odwrotności pierwiastka kwadratowego na każdej z kolejnych liczb - czyli wykonuję 209715200 takich operacji. Wynik każdej operacji zapisuję w innej tablicy niż tablica wejściowa. Testuję kolejno każdą z metod. Czas jest liczony w sekundach. Błąd jest liczony w porównaniu z metodą na FPU.
Wyniki na laptopie z Pentium 4-M 2GHz, DDR-266:
Input: 200x1048576 positive floats (max: 999.998291)
Normal time : 9.333797
Fast time : 2.626831 (speedup: 3.553254)
FastExtra time: 4.312687 (speedup: 2.164265)
SSE time : 1.755530 (speedup: 5.316797)
NORMAL/FAST:
mean error: -0.0000602813
error variance: 0.0000000182
max error: -0.0563874133
NORMAL/FASTEXTRA:
mean error: -0.0000001185
error variance: 0.0000000000
max error: -0.0001132605
NORMAL/SSE:
mean error: 0.0000000009
error variance: 0.0000000002
max error: -0.0032330328
Wyniki na desktopie z AMD Athlon X2 BE-2300, DDR2-800:
Input: 200x1048576 positive floats (max: 999.998291)
Normal time : 6.598592
Fast time : 1.681536 (speedup: 3.924146)
FastExtra time: 2.198866 (speedup: 3.000907)
SSE time : 0.680382 (speedup: 9.698363)
NORMAL/FAST:
mean error: -0.0000602813
error variance: 0.0000000182
max error: -0.0563874133
NORMAL/FASTEXTRA:
mean error: -0.0000001185
error variance: 0.0000000000
max error: -0.0001132605
NORMAL/SSE:
mean error: 0.0000014102
error variance: 0.0000000001
max error: -0.0032330328
Wersja SSE jest wyraźnie najszybsza, ale największy błąd wyniósł nawet 3e-3. Warto zwrócić uwagę na to jak różny jest średni błąd na różnych maszynach, ale to chyba normalne przy SSE - będzie to wyraźniejsze w drugim teście.
TEST: Normalize(float3)
Ponieważ odwrotność pierwiastka wykorzystuję przede wszystkim do normalizowania trójwymiarowych wektorów, postanowiłem jeszcze przetestować jak każda z tych metod się do tego nadaje. Tutaj obliczając błąd sprawdzam jak daleka jest długość wyjściowych wektorów od pożądanej jedynki.
Wyniki na Pentium 4-M 2GHz, DDR-266:
Input: 200x1048576 3d non zero vectors
Normal time : 10.012473
Fast time : 6.510379 (speedup: 1.537925)
FastExtra time: 7.751822 (speedup: 1.291628)
SSE time : 5.167708 (speedup: 1.937508)
NORMAL:
mean error: 0.0000000000
error variance: 0.0000000000
max error: -0.0000000503
FAST:
mean error: 0.0009001758
error variance: 0.0000003486
max error: 0.0017512076
FASTEXTRA:
mean error: 0.0000017375
error variance: 0.0000000000
max error: 0.0000046443
SSE:
mean error: -0.0000003470
error variance: 0.0000000133
max error: -0.0003259508
Wyniki na AMD Athlon X2 BE-2300, z DDR2-800:
Input: 200x1048576 3d non zero vectors
Normal time : 7.214732
Fast time : 4.266968 (speedup: 1.690833)
FastExtra time: 5.219265 (speedup: 1.382327)
SSE time : 2.117634 (speedup: 3.406978)
NORMAL:
mean error: 0.0000000000
error variance: 0.0000000000
max error: -0.0000000503
FAST:
mean error: 0.0009001758
error variance: 0.0000003486
max error: 0.0017512076
FASTEXTRA:
mean error: 0.0000017375
error variance: 0.0000000000
max error: 0.0000046443
SSE:
mean error: -0.0000209675
error variance: 0.0000000057
max error: -0.0002579081
Przyspieszenie w przypadku SSE jest zadowalające - na szybszej maszynie wręcz bardzo. Na laptopie (Pentium 4M) wersja Fast jest całkiem blisko SSE, chociaż znowu ma wyraźnie gorszą dokładność. Jednak tam zyski z SSE najprawdopodobniej giną przesłonięte przez opóźnienia wywołane fatalnie wolną pamięcią.
Jak to zwykle przy SSE, ważne jest odpowiednie ułożenie danych tak by dało się obliczać cztery wektory na raz, co w niektórych warunkach może być kłopotliwe, jeżeli nie niemożliwe. Dla wszystkich metod ułożyłem dane wektor po wektorze - tzn. xyzxyzxy..., za to dla SSE dane są wcześniej przekształcone w ciąg: xxxxyyyyzzzzxxx... . We wszystkich przypadkach pamięć jest wyrównana do granicy 16-bajtowej. Jeżeli idzie o implementację dla SSE, to wygląda następująco:
// SIZE - ilość wektorów
float xyz4[3*SIZE]; // wektory wejściowe (xxxxyyyyzzzzxxx...) (16-byte aligned)
float xyzs[3*SIZE]; // wektory wyjściowe - znormalizowane (16-byte aligned)
...
for (int i = 0; i < 3*SIZE; i += 12) {
__m128 r0 = _mm_load_ps(xyz4 + i);
__m128 r1 = _mm_load_ps(xyz4 + i + 4);
__m128 r2 = _mm_load_ps(xyz4 + i + 8);
__m128 r3 = r0;
__m128 r4 = r1;
__m128 r5 = r2;
r0 = _mm_mul_ps(r0, r0);
r1 = _mm_mul_ps(r1, r1);
r2 = _mm_mul_ps(r2, r2);
r0 = _mm_add_ps(r0, r1);
r0 = _mm_add_ps(r0, r2);
r0 = _mm_rsqrt_ps(r0);
r3 = _mm_mul_ps(r3, r0);
r4 = _mm_mul_ps(r4, r0);
r5 = _mm_mul_ps(r5, r0);
_mm_store_ps(xyzs + i, r3);
_mm_store_ps(xyzs + i + 4, r4);
_mm_store_ps(xyzs + i + 8, r5);
}
Mięso w pętli jest tłumaczone przez gcc na:
movl -24(%ebp), %eax
movaps (%eax,%ebx), %xmm3
leal 16(%ebx), %edx
movaps (%eax,%edx), %xmm4
leal 32(%ebx), %ecx
movaps (%eax,%ecx), %xmm5
movaps %xmm3, %xmm0
mulps %xmm3, %xmm0
movaps %xmm4, %xmm1
mulps %xmm4, %xmm1
movaps %xmm5, %xmm2
mulps %xmm5, %xmm2
addps %xmm1, %xmm0
addps %xmm2, %xmm0
rsqrtps %xmm0, %xmm0
mulps %xmm0, %xmm3
mulps %xmm0, %xmm4
mulps %xmm0, %xmm5
movl -40(%ebp), %eax
movaps %xmm3, (%eax,%ebx)
movl -40(%ebp), %eax
movaps %xmm4, (%eax,%edx)
movl -40(%ebp), %eax
movaps %xmm5, (%eax,%ecx)
Na koniec: 1.f/sqrtf(0.f)
Przy okazji normalizacji interesującym było jeszcze obejrzeć jak się te metody zachowują dla wektorów zerowych. FPU i SSE jako wynik 1/sqrt(0) zwracają inf (FPU bardzo powoli - można to jakoś przełączyć?), co dalej przemnożone przez 0.f daje nan. Za to w przypadku metod Fast i FastExtra dostałem bardzo duże wartości (Fast: 1.98e+19, FastExtra: 2.97e+19), które śmiało można przemnożyć przez 0. Wynik będzie taki, że normalizacja wektora na FPU i SSE daje wektor [NaN, NaN, NaN], za to pozostałe przybliżone metody dadzą [0.f, 0.f, 0.f]. Zależy od zastosowania, który wynik jest bardziej pożądany, chociaż pewno w większości przypadków lepiej jest standardowo dorzucić ifa po obliczeniu sumy kwadratów współrzędnych, ale może to skomplikować obliczenia na SSE.
Błądy ujemny i dodatnie:
Sprawdziłem jeszcze (na Pentium 4M), dla normalizacji wektora, jak dokładnie wyglądają najgorsze błędy dodatnie i ujemne:
NORMAL:
max +error: 0.0000000502
max -error: -0.0000000503
FAST:
max +error: 0.0017512076
max -error: -0.0000000384
FASTEXTRA:
max +error: 0.0000046443
max -error: -0.0000000468
SSE:
max +error: 0.0003250295
max -error: -0.0003259508
Jak widać błąd ujemny jest minimalny - nawet mniejszy niż przy FPU. Błąd ujemny wskazuje na to o ile długość znormalizowanego wektora przekracza jedynkę. Czyli w przypadku metod Fast i FastExtra większość błędów powoduje wyprodukowanie wektorów trochę krótszych niż wektory normalne. Przy tak dużym błędzie absolutnym jest to fantastyczna własność - przynajmniej w moich zastosowaniach.
Starsze posty dostępne w archiwum
Ładuję...
