szybkie 1/sqrt, sqrt i dzielenie

Obliczenie 1/sqrt(a):
najpierw sqrt i jeszcze dzielenie - obie operacje są długo wykonywane na FPU, i razem to trwa ze 100 taktów.

Za to mnożenie i dodawanie są bardzo szybkie, więc może da się szybciej to wyliczyć.
Wyrażenia typu: q / sqrt(x^{2 + y}2) w wielu problemach bardzo często obliczamy, więc gdyby to przyspieszyć z 5 razy byłoby znacznie lepiej.

Tu jest taki zabawny algorytm:
http://en.wikipedia.org/wiki/Fast_inverse_square_root

Trochę mała precyzja, chyba z 0.2% błędu (z jedną iteracją), czyli z 9 bitów.
Żeby uzyskać pełną dokładność należałoby wykonać kolejne iteracje Newtona - chyba minimum dwie, dla precyzji float, a dla double jeszcze jedną albo i dwie.

Powiedzmy, że mamy dobre 8 bitów, a w każdym kroku Newtona to się praktycznie podwaja:
-> 16 -> 32 -> 64; tylko 3 dodatkowe iteracje dla 64 bitów?

Razem są 4 iteracje typu:
x = x*(1.5 - a2xx); gdzie: a2 = a*0.5;

Są tu 3 mnożenia i jedno odejmowanie, razem 4.
4 x 4 = 16 operacji, co pewnie byłoby szybsze od tej bezpośredniej metody, ale raczej niewiele.

Trzeba zmniejszyć liczbę iteracji, przynajmniej o jedną, a najlepiej o dwie.

Nie można przyspieszyć zbieżność metody Newtona?
Wiemy jak błąd tu maleje, więc wyliczamy go (jedna iteracja) i coś tam korygujemy przed kolejnym krokiem... ale jak?

100 taktów ale na jakim procesorze? Wszystkie instrukcje, które cię interesują znajdują się w zestawie instrukcji SSE dla liczb pojedynczej precyzji, instrukcje dla podwójnej prezycji są w SSE2.

Dla przykładu na Intel Sandy Bridge FSQRT dla long double (80-bit) to max 24 takty, tak samo FDIV.

A co jest nie tak z SSE / rsqrtss?

Przykład:
http://social.msdn.microsoft.com/Forums/en-US/vclanguage/thread/99fa1e9f-9dc3-46a0-9c36-088762d4f1c9/

Z SSE pewnie byłoby najszybciej w przypadku 1/sqrt.
Ale rozkaz rsqrt oblicza tylko przybliżenie - precyzja 12 bitów, czyli to jest coś w stylu tego invSqrt z Quake.

Poza tym, nie należy zakładać tu na beton precyzji double,
lecz przynajmniej quad; w tym przypadku SSE da nam tylko przybliżenie.

No to zrób sobie invSqrt na FPU, a potem skonwertuj na quada i dalej przybliżaj.

Raczej niewiele odkryłeś.

Sandy Bridge ma 24 takty na fsqrt i fdiv, ale fmul i fadd proporcjonalnie mniejsze - pewnie 1 z optymalnym kodem.

Dalej.
sqrt: mając 1/sqrt(a) = a^{-0.5, mnożymy tylko: a * a}-0.5 = sqrt(a) i gotowe.

Dzielenie. Metoda Newtona dla równania: 1/x - a, daje taki schemat:
x = x*(2 - a*x);

Wystarczy przybliżenie 8 bitów, a wtedy 3 iteracje i mamy 64 bity.
3 operacje na iterację, więc tu też powinno być szybciej... w procesorach pewnie tak samo to liczą, więc dziwne, że to tak długo trwa).

Sprawdźmy jak to działa: a = 1 + m; i m z przedziału <0, 1);
zatem szukane: 1/a = 1/(1+m) = 1 - m + m^{2 - m}3 + ...

Możemy startować od: x = 1 - m, 50% błędu dla m = 0.5;
a = 1.5; 1/a = 2/3 = 0.6666...

x = (1-m)(2 - (1+m)(1-m)) = (1-m)(2 - (1-m^{2)) = (1-m)(1 + m}2) = 5/8 = 0.625;
x = 5/8 (2 - 3/2 * 5/8) = 0.6640625;
x = 0.666656494140625

a gdy wystartujemy z: x = 0.66, czyli 1% błędy:
x = 0.66 * (2 - 1.5*.66) = 0.6666
x = 0.66666666;

no, ale przydałoby się przyspieszenie zbieżności.
.........

Ale przecież: (a^-0.5)2 = 1/a;
Nie potrzeba bawić się z dzieleniem - jeden rsqrt załatwia tu wszystko.

Agner Fog podaje (dla Intel Sandy Bridge):

Jak widać jest kilka jednostek do FMULowania i przy dobrym ustawieniu instrukcji można je przetworzyć potokowo uzyskując ten jeden cykl na instrukcję, ale w przypadku metody Newtona ciężko to tak spotokować, bo przecież kolejne wyniki zależą od siebie i nie można ich oddzielnie przetworzyć.

Przy SSE na Sandy Bridge dodawanie i mnożenie ma takie same Latency i Reciprocal Throughput, natomiast dzielenie, wyciąganie pierwiastka czy odwrotności jest szybsze, ponieważ mamy tutaj do czynienia z niższą precyzją (32-bit float lub 64-bit double zamiast 80-bit long double).

AMD używa dla przykładu metody Goldschmidta: http://en.wikipedia.org/wiki/Division_(digital)#Goldschmidt_division
Wszelkie próby uzyskania prędkości wyższej niż używając bezpośrednio instrukcji do dzielenia, wyciągania pierwiastka czy odwrotności tak aby uzyskać taką samą dokładność są wg mnie z góry skazane na niepowodzenie, no chyba że klecisz coś na jakieś bardzo słabe lub stare procesory. Dzielenie to bardzo typowa operacja i była ona optymalizowana w prockach od dekad, myślenie że uda się oklepanym od stuleci Newtonem osiągnąć coś lepszego to totalna pomyłka.

@Wii: Pytanie czego od nas oczekujesz?
Implementacji swojego algorytmu, czy podanie gotowych rozwiązań?

Bo "szybkich" teoretycznych algorytmów jest kilka:
http://www.isrn.com/journals/sp/2011/615934/
http://bit.ly/qe6Sv5
http://bit.ly/pQTOJC

A może zamiast szukać wydajności pojedynczego pierwiastka przemyślisz algorytm i np. zastosujesz wektory?
http://intel.ly/qVxUBR

Wibowit napisał(a)

AMD używa dla przykładu metody Goldschmidta: http://en.wikipedia.org/wiki/Division_(digital)#Goldschmidt_division
Wszelkie próby uzyskania prędkości wyższej niż używając bezpośrednio instrukcji do dzielenia, wyciągania pierwiastka czy odwrotności tak aby uzyskać taką samą dokładność są wg mnie z góry skazane na niepowodzenie, no chyba że klecisz coś na jakieś bardzo słabe lub stare procesory. Dzielenie to bardzo typowa operacja i była ona optymalizowana w prockach od dekad, myślenie że uda się oklepanym od stuleci Newtonem osiągnąć coś lepszego to totalna pomyłka.

Sam sobie zaprzeczasz: starsze procesory potrzebowały z 40-80 taktów na dzielenie i fsqrt, a dopiero te najnowsze 24, więc jakoś marnie pracowali nad tymi instrukcjami, a instrukcji 1/sqrt (obliczającej dokładnie) w ogóle nie ma, więc tu 2x tracisz.

Można było zwiększyć szybkość już na Pentium, czyli z 10 lat temu.

http://www.intel-assembler.it/portale/5/fpu-timing/8087-pentium-coprocessor-timing-pairing.asp

Pentium: fdiv 40, fsqrt 80, a mnoży i dodaje wielokrotnie szybciej.
dla 1/sqrt otrzymasz tu 120.

Ten algorytm Goldsmidta nie jest raczej szybszy od Newtona.
Wystarczy odwrotność części ułamkowej: 0 < x < 1;
co można zapisać: x = 1 - m;

i tam jest liczone coś takiego: (1 + m)(1 + m²⁾(1 + m</sup>4)(1 + m^8) ...
zbieżność kwadratowa, czyli tak samo jak Newton.

x = .5; obliczamy: 1/x = 2;
(1+1/2) = 1.5
1.5*(1 + 1/4) = 1.875
1.875*(1+1/16) = 1.9921875
1.9921875*(1 + 1/256) = 1.999969482421875

Ty chcesz optymalizować dla staroci czy dla najnowszych procesorów? Dla staroci być może się opłaca, dla nowszych już nie.

Możliwe, że AMD wybrało Goldschmidta dlatego, że daje się zaimplementować z użyciem mniejszej liczby tranzystorów. Innym powodem może być to, że mnożenia w Goldschmidcie lepiej się zrównoleglają, np mając policzone (1 + x^{2) oraz x}2 można jednocześnie dokonać dwóch mnożeń: (1 + x) * (1 + x^{2) oraz x}2 * x^{2 = x}4, potem jednocześnie kolejne dwa mnożenia: (wynik z poprzedniego) * (1 + x^{4) oraz x}4 * x^{4 = x}8.

Dzielenie na FPU było tak czy siak dość szybkie nawet już na Pentium Plain, tam było <40 taktów na dzielenie w pełnej precyzji (tzn 80-bit long double). Optymalizacje polegają generalnie na stworzeniu bardziej skomplikowanych, ale wymagających mniej kroków algorytmów, użyciu większych lookup tables, etc Zmniejsza to budżet tranzystorów na inne funkcjonalności, więc inżynierowie starają się mocno optymalizować tylko te najbardziej potrzebne instrukcje.

Możesz sobie np stablicować pierwiastki dla np pierwszych 20 bitów mantysy i wtedy pewnie dostaniesz coś szybszego niż to co jest wbudowane w procka, no ale z drugiej strony obciążysz znacznie pamięć podręczną + kontroler pamięci, jeżeli będziesz chciał pierwiastkować różne liczby.

Gdyby się dało prosto przyspieszyć dzielenie bez drastycznego zwiększania liczby potrzebnych tranzystorów to inżynierowie Intela dawno by to zrobili, nie są w ciemię bici.

Zamiast przyspieszać dzielenie i pierwiastkowanie oddzielnie, zrób szybką wersję instrukcji, która robi obydwie rzeczy naraz.

Wibowit napisał(a)

Ty chcesz optymalizować dla staroci czy dla najnowszych procesorów? Dla staroci być może się opłaca, dla nowszych już nie.

Odwrotnie.
Na starociach nie można przyspieszyć, bo tam mnożenie i dodawanie jest też wolne.
Tylko te nowsze można przyspieszać, i pewnie łącznie z tymi, które wyprodukują w 2012r.

Wibowit napisał(a)

Możesz sobie np stablicować pierwiastki dla np pierwszych 20 bitów mantysy i wtedy pewnie dostaniesz coś szybszego niż to co jest wbudowane w procka, no ale z drugiej strony obciążysz znacznie pamięć podręczną + kontroler pamięci, jeżeli będziesz chciał pierwiastkować różne liczby.
Gdyby się dało prosto przyspieszyć dzielenie bez drastycznego zwiększania liczby potrzebnych tranzystorów to inżynierowie Intela dawno by to zrobili, nie są w ciemię bici.

Pewnie ci fachowcy planują tablicować 20 bitów, ale pamięci zbyt drogie, więc czekają.

Odwrotnie.
Na starociach nie można przyspieszyć, bo tam mnożenie i dodawanie jest też wolne.
Tylko te nowsze można przyspieszać, i pewnie łącznie z tymi, które wyprodukują w 2012r.

No to spróbuj waćpan, albo przynajmniej poczytaj timingi instrukcji.

Na PPlain:

Im nowszy procesor tym trudniej jest przyspieszyć dzielenie w przypadku ogólnym i tym mniejszy jest sens czynienia tego.

Not much separates Penryn from the first generation of Core 2 architecture-wise, but Intel has taken the opportunity to make a few improvements. We present some of the most notable changes below.
Fast Radix-16 divider and Super Shuffle Engine
Fast Radix-16 divider is an improvement of the general division instruction. Earlier techniques worked on two bits at a time while the processor now has the possibility to work on 4 bits at a time, which results in twice the performance when dividing. A prominent advantage of changes to basic instructions is that these affect already existing software which uses division, as opposed to special instructions which require a recompile of the program with an updated compiler, or be manually updated.

Wibowit napisał(a)

No to spróbuj waćpan, albo przynajmniej poczytaj timingi instrukcji.

http://www.friweb.hu/instlatx64/
Przeczytałem.

Wibowit napisał(a)

Na PPlain:
||=Instruction||Latency||Reciprocal Throughput
||FADD||3||1
||FMUL||3||ok 1.5 (zależy czy kolejna instrukcja to też FMUL)
||FDIV||39 dla przecyzji 80-bitów (long double)||-

No i tak jest na pentium 4 i dalej. Dopiero teraz zmniejszyli dwa razy (powiększyli lookup table, bo pamięci tańsze).

Wibowit napisał(a)

Im nowszy procesor tym trudniej jest przyspieszyć dzielenie w przypadku ogólnym i tym mniejszy jest sens czynienia tego.

Nie ma sensu, ale dla intela i amd.
Taktyka znana - podobnie Sergiej Bubka robił.

FPU i SSE działają niezależnie (równolegle)?

Zamiast wymyślań bajki o końcu świata (już wszystko zrobiono, więc co pozostało?), lepiej podaj metodę przyspieszenia zbieżności.

Np. Halley: x = x/8 * [15 - q(10-3q)]; q = ax^2; metoda 3-go rzędu dla 1/sqrt(a).

log_a(x)/log_b(x) = log_a(b), a więc jeśli przechodzisz z rzędu zbieżności a do rzędu zbieżności b, ale wykonujesz w każdym kroku co najmniej log_a(b) razy więcej pracy to otrzymujesz bardziej pracochłonny algorytm.

FPU i SSE działają niezależnie (równolegle)?

Zależy od procesora. Nawet jeśli na danym procesorze działają równolegle, to możesz być ograniczony przez dekoder instrukcji, dekodowanie instrukcji x86 jest powolne i prądożerne.

Intel używa w swoich prockach czegoś na wzór dzielenia szkolnego, tyle że nieco mniej operacji na krok (non-restoring division zamiast restoring division) oraz używa lookup tables. Przy Radix-16 w nowych procesorach Intela liczone są 4 bity w jednym kroku, a więc na mantysę 64-bit (w 80-bit long double) jest 16 kroków, dla double 14 kroków, dla float 6 kroków, gdzie jeden krok jest robiony w jednym takcie.

W algorytmie Halleya na jeden krok jest dość dużo operacji, zwłaszcza, że operacja mnożenia jest dużo kosztowniejsza do zaimplementowania w hardware niż operacja dodawania.

Metody Halleya czy Newtona są całkowicie "seryjne" (serial, tzn instrukcje muszą być wykonywane in-order) natomiast przy Goldschmidcie można trochę wykorzystać superskalarność (tzn można wykonywać instrukcje out-of-order).
Jeśli jeden krok jest dwa razy szybszy, to tak jakbyś miał dwa razy mniej kroków :p

Wibowit napisał(a)

log_a(x)/log_b(x) = log_a(b), a więc jeśli przechodzisz z rzędu zbieżności a do rzędu zbieżności b, ale wykonujesz w każdym kroku co najmniej log_a(b) razy więcej pracy to otrzymujesz bardziej pracochłonny algorytm.

Nieprawda.
Patrz na dzielenie z Newtona:

x = x + x(1 - ax);
składasz to dwa razy ze sobą i będzie metoda rzędu 4, ale najpierw robimy tak:
p = 1 - ax; i tylko raz to obliczamy:
zatem: x = x(1 + p);

dopiero teraz składamy i otrzymasz: x = x(1 + p)(1 + p^2);
jeszcze raz i będzie 8 rzędu: x = x(1 + p)(1 + p^{2)(1 + p}4);

I to jest chyba algorytm tego Golfshmidta (albo coś podobnego).

To jest zwyczajny szereg Taylora:
1/a = 1/(1/x - (a - 1/x)) = x/(1 - p) = x(1 + p + p^{2 + ... ) = x(1+p)(1+p}2)...
x - to przybliżenie początkowe 1/a; wystarczy 8 bitów, czyli tablica 256 bajtów.

rsqrt podobnie, ale tu nie faktoryzuje się tak łatwo ten szereg:
1/sqrt(a) = x(1 + 1/2 p + 3/8 p^{2 + 8/16 p}3 + ...); p = ...

Wibowit napisał(a)

Metody Halleya czy Newtona są całkowicie "seryjne" (serial, tzn instrukcje muszą być wykonywane in-order) natomiast przy Goldschmidcie można trochę wykorzystać superskalarność (tzn można wykonywać instrukcje out-of-order).
Jeśli jeden krok jest dwa razy szybszy, to tak jakbyś miał dwa razy mniej kroków :p

Masz lepiej niż Gold...tegos.
Nawet programując na CPU będzie szybciej.
A używając na starcie przybliżenia 12 bit z SSE: 12 * 5 = 60, metoda rzędu 4.5 wystarczy.