Dlaczego skalarny parametr sqrt (x) SSE jest wolniejszy niż rsqrt (x) * x?

Profilowałem część naszej podstawowej matematyki na Intel Core Duo i patrząc na różne podejścia do pierwiastka kwadratowego zauważyłem coś dziwnego: używając operacji skalarnych SSE, szybciej jest wziąć odwrotność pierwiastka kwadratowego i pomnożyć go aby uzyskać sqrt, niż użyć natywnego kodu operacji sqrt!

Testuję to z pętlą coś takiego:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

Próbowałem tego z kilkoma różnymi ciałami dla TestSqrtFunction i mam kilka czasów, które naprawdę drapią mnie po głowie. Zdecydowanie najgorsze było użycie natywnej funkcji sqrt () i pozwolenie „inteligentnemu” kompilatorowi na „optymalizację”. Przy 24ns / float, przy użyciu FPU x87 było to żałośnie złe:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

Następną rzeczą, jaką próbowałem, było użycie funkcji wewnętrznej, aby zmusić kompilator do użycia skalarnego kodu operacji sqrt SSE:

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

To było lepsze, przy 11,9ns / float. Wypróbowałem również zwariowaną technikę przybliżenia Newtona-Raphsona Carmacka , która działała nawet lepiej niż sprzęt, przy 4,3 ns / float, chociaż z błędem 1 na 2 ¹⁰ (co jest zbyt duże dla moich celów).

Doozy miał miejsce, gdy próbowałem opcją SSE dla odwrotności pierwiastka kwadratowego, a następnie użyłem mnożenia, aby uzyskać pierwiastek kwadratowy (x * 1 / √x = √x). Mimo to trwa dwie operacje zależne było najszybsze rozwiązanie zdecydowanie na 1.24ns / pływaka i dokładnością do 2 ^-14 :

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Moje pytanie brzmi: co daje ? Dlaczego kod operacji pierwiastka kwadratowego wbudowany w sprzęt SSE jest wolniejszy niż jego synteza z dwóch innych operacji matematycznych?

Jestem pewien, że to tak naprawdę koszt samej operacji, bo zweryfikowałem:

• Wszystkie dane mieszczą się w pamięci podręcznej, a dostęp jest sekwencyjny
• funkcje są wbudowane
• rozwijanie pętli nie robi różnicy
• flagi kompilatora są ustawione na pełną optymalizację (a montaż jest dobry, sprawdziłem)

( edytuj : stephentyrone poprawnie wskazuje, że operacje na długich ciągach liczb powinny wykorzystywać wektoryzację operacji spakowanych w SIMD, na przykład rsqrtps- ale struktura danych tablicy jest tutaj tylko do celów testowych: to, co naprawdę próbuję zmierzyć, to wydajność skalarna do użycia w kodzie których nie można wektoryzować).

sqrtssdaje poprawnie zaokrąglony wynik. rsqrtsspodaje przybliżenie odwrotności, z dokładnością do około 11 bitów.

sqrtssgeneruje znacznie dokładniejsze wyniki, gdy wymagana jest dokładność. rsqrtssistnieje w przypadkach, gdy wystarczy przybliżenie, ale wymagana jest prędkość. Jeśli przeczytasz dokumentację Intela, znajdziesz również sekwencję instrukcji (odwrotne przybliżenie pierwiastka kwadratowego, po którym następuje pojedynczy krok Newtona-Raphsona), która zapewnia prawie pełną precyzję (~ 23 bity dokładności, jeśli dobrze pamiętam) i nadal jest nieco szybciej niż sqrtss.

edycja: Jeśli szybkość ma kluczowe znaczenie i naprawdę wywołujesz to w pętli dla wielu wartości, powinieneś używać wektoryzowanych wersji tych instrukcji rsqrtpslub sqrtpsobie przetwarzają cztery zmiennoprzecinkowe na instrukcję.

Dotyczy to również podziału. MULSS (a, RCPSS (b)) jest znacznie szybszy niż DIVSS (a, b). W rzeczywistości jest nadal szybszy, nawet jeśli zwiększysz jego precyzję za pomocą iteracji Newtona-Raphsona.

Intel i AMD zalecają tę technikę w swoich podręcznikach optymalizacji. W aplikacjach, które nie wymagają zgodności ze standardem IEEE-754, jedynym powodem używania div / sqrt jest czytelność kodu.

Zamiast udzielać odpowiedzi, która w rzeczywistości może być niepoprawna (nie zamierzam też sprawdzać ani dyskutować o pamięci podręcznej i innych rzeczach, powiedzmy, że są identyczne) spróbuję wskazać źródło, które może odpowiedzieć na twoje pytanie.
Różnica może polegać na sposobie obliczania sqrt i rsqrt. Więcej informacji można znaleźć tutaj http://www.intel.com/products/processor/manuals/ . Proponuję zacząć od przeczytania o funkcjach procesora, których używasz, jest trochę informacji, szczególnie o rsqrt (procesor używa wewnętrznej tabeli przeglądowej z ogromnym przybliżeniem, co znacznie ułatwia uzyskanie wyniku). Może się wydawać, że rsqrt jest o wiele szybszy niż sqrt, że 1 dodatkowa operacja mul (co nie jest zbyt kosztowna) może nie zmienić sytuacji tutaj.

Edycja: Kilka faktów, o których warto wspomnieć:
1. Kiedyś robiłem mikro optymalizacje dla mojej biblioteki graficznej i użyłem rsqrt do obliczania długości wektorów. (zamiast sqrt pomnożyłem sumę kwadratu przez rsqrt, co jest dokładnie tym, co zrobiłeś w swoich testach) i wypadło lepiej.
2. Obliczenie rsqrt przy użyciu prostej tabeli przeglądowej może być łatwiejsze, ponieważ dla rsqrt, gdy x dochodzi do nieskończoności, 1 / sqrt (x) idzie do 0, więc dla małych x wartości funkcji się nie zmieniają (dużo), podczas gdy dla sqrt - dąży do nieskończoności, więc to taki prosty przypadek;).

Ponadto wyjaśnienie: nie jestem pewien, gdzie znalazłem to w książkach, do których linkowałem, ale jestem prawie pewien, że czytałem, że rsqrt używa jakiejś tabeli odnośników i powinno być używane tylko wtedy, gdy wynik nie musi być dokładne, chociaż - ja też mogę się mylić, tak jak to było jakiś czas temu :).

Newton-Raphson zbiega się do zera f(x)przy użyciu przyrostów równych -f/f' gdzie f'jest pochodną.

Ponieważ x=sqrt(y)możesz spróbować rozwiązać f(x) = 0za xpomocą f(x) = x^2 - y;

Wtedy przyrost jest następujący: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x który ma powolny podział.

Możesz wypróbować inne funkcje (np. f(x) = 1/y - 1/x^2), Ale będą one równie skomplikowane.

Spójrzmy 1/sqrt(y)teraz. Możesz spróbować f(x) = x^2 - 1/y, ale będzie to równie skomplikowane: dx = 2xy / (y*x^2 - 1)na przykład. Jednym z nieoczywistych alternatywnych opcji f(x)jest:f(x) = y - 1/x^2

Następnie: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

Ach! To nie jest trywialne wyrażenie, ale masz w nim tylko mnożenia, bez dzielenia. => Szybciej!

I: pełny krok aktualizacji new_x = x + dxbrzmi:

x *= 3/2 - y/2 * x * x co też jest łatwe.

Istnieje wiele innych odpowiedzi na to pytanie sprzed kilku lat. Oto, co jest słuszne w konsensusie:

• Instrukcje rsqrt * obliczają przybliżenie do odwrotności pierwiastka kwadratowego, dobre do około 11-12 bitów.
• Jest zaimplementowany z tablicą przeglądową (tj. ROM) indeksowaną przez mantysę. (W rzeczywistości jest to skompresowana tablica przeglądowa, podobna do starych tablic matematycznych, wykorzystująca korekty mniej znaczących bitów, aby zaoszczędzić na tranzystorach).
• Powodem, dla którego jest dostępny, jest to, że jest to wstępne oszacowanie używane przez FPU dla "prawdziwego" algorytmu pierwiastka kwadratowego.
• Istnieje również przybliżona wzajemna instrukcja, rcp. Obie te instrukcje są wskazówką, jak FPU implementuje pierwiastek kwadratowy i dzielenie.

Oto dlaczego konsensus się nie zgadzał:

• Jednostki FPU ery SSE nie używają Newtona-Raphsona do obliczania pierwiastków kwadratowych. To świetna metoda w oprogramowaniu, ale byłoby błędem wdrażanie jej w ten sposób w sprzęcie.

Algorytm NR do obliczania odwrotności pierwiastka kwadratowego ma ten krok aktualizacji, jak zauważyli inni:

x' = 0.5 * x * (3 - n*x*x);

To dużo mnożenia zależnego od danych i jedno odejmowanie.

Poniżej znajduje się algorytm, którego faktycznie używają nowoczesne jednostki FPU.

Biorąc pod uwagę b[0] = n, załóżmy, że możemy znaleźć szereg liczb Y[i], które b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2zbliżają się do 1. Następnie rozważ:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

Jasne x[n]podejście sqrt(n)i y[n]podejście 1/sqrt(n).

Możemy użyć kroku aktualizacji Newtona-Raphsona do odwrotności pierwiastka kwadratowego, aby uzyskać dobry Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Następnie:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

i:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

Następna kluczowa obserwacja jest taka b[i] = x[i-1] * y[i-1]. Więc:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Następnie:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

Oznacza to, że mając początkowe x i y, możemy użyć następującego kroku aktualizacji:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

Lub, nawet bardziej wyszukane, możemy ustawić h = 0.5 * y. To jest inicjalizacja:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

A to jest krok aktualizacji:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

To jest algorytm Goldschmidta, który ma ogromną zaletę, jeśli implementujesz go w sprzęcie: „pętla wewnętrzna” to trzy wielokrotne dodawanie i nic więcej, a dwa z nich są niezależne i można je potokować.

W 1999 r. Jednostki FPU potrzebowały już potokowego obwodu dodawania / odejmowania i potokowego obwodu wielokrotnego, w przeciwnym razie SSE nie byłoby zbyt „strumieniowe”. W 1999 r. Potrzebny był tylko jeden z każdego obwodu, aby zaimplementować tę wewnętrzną pętlę w sposób w pełni potokowy bez marnowania dużej ilości sprzętu tylko na pierwiastek kwadratowy.

Dziś oczywiście połączyliśmy mnożenie i dodawanie ujawnione programiście. Ponownie, pętla wewnętrzna to trzy potokowe FMA, które są (znowu) ogólnie przydatne, nawet jeśli nie obliczasz pierwiastków kwadratowych.

Jest to szybsze, ponieważ te instrukcje ignorują tryby zaokrąglania i nie obsługują wyjątków zmiennoprzecinkowych ani zdernormalizowanych liczb. Z tych powodów znacznie łatwiej jest potokować, spekulować i wykonać inne instrukcje fp.