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

Dotyczy to również podziału. MULSS(a,RCPSS(b)) jest o wiele szybszy niż DIVSS (a,b). W rzeczywistości jest jeszcze szybszy, nawet jeśli zwiększysz jego precyzję iteracją 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życia div / sqrt jest czytelność kodu.

Zamiast podawać odpowiedź, która faktycznie może być niepoprawna (Nie będę też sprawdzać ani dyskutować o cache i innych rzeczach, powiedzmy, że są identyczne) postaram się wskazać źródło, które może odpowiedzieć na twoje pytanie.
Różnica może leżeć w tym, jak sqrt i rsqrt są obliczane. Możesz przeczytać więcej tutaj http://www.intel.com/products/processor/manuals / . proponuję zacząć od czytania o funkcjach procesora, których używasz, są pewne informacje, zwłaszcza o rsqrt (cpu używa wewnętrznej tabeli wyszukiwania z dużym przybliżeniem, co znacznie ułatwia uzyskanie wyniku). Może się wydawać, że rsqrt jest o tyle szybszy od sqrt, że 1 dodatkowa operacja wielostopniowa (która nie jest zbyt kosztowna) może nie zmienić sytuacji.

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

Również wyjaśnienie: nie jestem pewien, gdzie znalazłem to w książkach, które linkowałem, ale jestem prawie pewien, że czytałem, że rsqrt jest używając jakiejś tabeli lookup, i powinna być używana tylko wtedy, gdy wynik nie musi być dokładny , chociaż-mogę się mylić, jak to było jakiś czas temu :).

Newton-Raphson jest zbieżny do zera z f(x) używając przyrostów równych -f/f' Gdzie f' jest pochodną.

Dla x=sqrt(y), Możesz spróbować rozwiązać f(x) = 0 dla x używając f(x) = x^2 - y;

Wtedy przyrost wynosi: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x która ma w sobie powolny podział.

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

Spójrzmy na 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. Jedna nieoczywista alternatywa wybór dla f(x) to: f(x) = y - 1/x^2

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

I: pełny krok aktualizacji new_x = x + dx następnie brzmi:

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

Istnieje wiele innych odpowiedzi na to już od kilku lat. Oto co konsensus ma rację:

• instrukcje rsqrt * obliczają przybliżenie do wzajemnego pierwiastka kwadratowego, dobre do około 11-12 bitów.
• jest zaimplementowany z tabelą wyszukiwania (np. ROM) indeksowaną przez mantysę. (W rzeczywistości jest to skompresowana tabela wyszukiwania, podobna do starych tabel matematycznych, wykorzystująca korekty bitów niskiego rzędu, aby zapisać na tranzystorach.)
• The powodem, dla którego jest on dostępny, jest to, że jest to wstępne oszacowanie używane przez FPU dla" rzeczywistego " algorytmu pierwiastka kwadratowego.
• jest też przybliżona Instrukcja wzajemna, rcp. Obie te instrukcje są wskazówką jak FPU implementuje pierwiastek kwadratowy i podział.

Oto co konsensus się pomylił:

• SSE-era FPU nie używa Newtona-Raphsona do obliczania pierwiastków kwadratowych. Jest to świetna metoda w oprogramowaniu, ale błędem byłoby zaimplementować ją w ten sposób w sprzęt.

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 współczesne FPU faktycznie używają.

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

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

Wyraźnie x[n] podejść sqrt(n) i y[n] podejść 1/sqrt(n).

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

Potem:

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

I:

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

Kolejną kluczową obserwacją jest to, że 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])

Potem:

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

To znaczy, biorąc pod uwagę 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 fantazyjne, możemy ustawić h = 0.5 * y. Jest to inicjalizacja:

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

A to jest aktualizacja krok:

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

Jest to algorytm Goldschmidta, który ma ogromną zaletę, jeśli implementujesz go w sprzęcie: "wewnętrzna pętla" to trzy mnożniki i nic więcej, a dwa z nich są niezależne i mogą być pipelinowane.

W 1999 r.FPU już potrzebowało pipelined obwodu dodawania/odejmowania i Pipelined obwodu mnożenia, w przeciwnym razie SSE nie byłoby zbyt "strumieniowe". Tylko jeden z każdego obwodu był potrzebny w 1999 roku, aby wdrożyć tę wewnętrzną pętlę w pełni pipelined sposób bez marnowanie dużo sprzętu tylko na pierwiastek kwadratowy.

Dzisiaj, oczywiście, połączyliśmy multiply-add z programistą. Ponownie, wewnętrzna pętla to trzy pipelinowe FMA, które są (ponownie) ogólnie użyteczne, 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 punktów pływających lub liczb dernormalizowanych. Z tych powodów o wiele łatwiej jest rurociąg, spekulować i wykonywać inne instrukcje fp Nie w porządku.