Lowlevel

Lowlevel => Softwareentwicklung => Thema gestartet von: TomCat am 16. February 2017, 15:10

Titel: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 16. February 2017, 15:10
hallo,
ich möchte die Quadratwurzel aus einer 32-bit integer Zahle berechnen.
Z.B. aus 50 wäre die Wurzel dann 7.

Mir geht es hier nicht um Genauigkeit. Eine Abweichung von 10-15 % vom eigentlichen Wert nehme ich in Kauf.
Es geht nur um maximale Geschwindigkeit !!
Wie würde man da am besten vorgehen?

THX
TomCat
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: kevin am 16. February 2017, 17:58
Naja, normal nimmt man da wohl irgendein iteratives Verfahren, z.B. das Newtonverfahren. Wenn du schnell und wenig präzise willst, machst du halt entsprechend weniger Iterationen.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 17. February 2017, 10:52
naja das iterationsverfahren muss schon 3 bis 4 mal durchgeführt werden. Das sind einige Maschinenbefehle die da durchgejagt werden müssen. Blöderweise auch Divisionen die relativ langsam sind. Das wird dann langsamer als gleich die FPU zu bemühen.

 :-(
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: kevin am 17. February 2017, 12:21
Wenn du einen Float kriegen kannst (und dann auf der Binärdarstellung davon mit Integeroperationen arbeiten), ist eventuell das hier (https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Approximations_that_depend_on_the_floating_point_representation) was für dich. Soll dann am Ende irgendwie so aussehen:
val_int = (1 << 29) + (val_int >> 1) - (1 << 22) + a;
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: Jidder am 18. February 2017, 21:55
Ich hätte da auch noch einen Vorschlag (für Fehler < 10%):

if (x >= 2875143395) return 58982;
if (x >= 1924682769) return 48258;
if (x >= 1288424002) return 39484;
if (x >= 862498712) return 32305;
if (x >= 577375171) return 26431;
if (x >= 386507345) return 21625;
if (x >= 258736322) return 17693;
if (x >= 173203653) return 14476;
if (x >= 115946247) return 11844;
if (x >= 77616909) return 9691;
if (x >= 51958427) return 7929;
if (x >= 34782087) return 6487;
if (x >= 23283876) return 5307;
if (x >= 15586727) return 4342;
if (x >= 10434090) return 3553;
if (x >= 6984804) return 2907;
if (x >= 4675778) return 2378;
if (x >= 3130066) return 1946;
if (x >= 2095333) return 1592;
if (x >= 1402661) return 1302;
if (x >= 938971) return 1065;
if (x >= 628567) return 872;
if (x >= 420776) return 713;
if (x >= 281677) return 583;
if (x >= 188560) return 477;
if (x >= 126226) return 390;
if (x >= 84498) return 319;
if (x >= 56565) return 261;
if (x >= 37866) return 214;
if (x >= 25348) return 175;
if (x >= 16968) return 143;
if (x >= 11359) return 117;
if (x >= 7604) return 95;
if (x >= 5090) return 78;
if (x >= 3407) return 64;
if (x >= 2281) return 52;
if (x >= 1527) return 42;
if (x >= 1022) return 35;
if (x >= 684) return 28;
if (x >= 458) return 23;
if (x >= 306) return 19;
if (x >= 205) return 15;
if (x >= 137) return 12;
if (x >= 91) return 10;
if (x >= 61) return 8;
if (x >= 41) return 7;
if (x >= 27) return 5;
if (x >= 18) return 4;
if (x >= 12) return 3;
if (x >= 8) return 3;
if (x >= 5) return 2;
if (x >= 3) return 2;
if (x >= 2) return 1;
if (x >= 1) return 1;
if (x >= 1) return 1;
if (x >= 0) return 0;

Der Fehler für x unter 1000 ist zwar teilweise größer als 10% und da sind redundante if-Abfragen drin. Das liegt vermutlich am Algorithmus, mit dem ich den Code generiert habe. Allerdings sind das auch weniger als 0,00001% des Definitionsbereichs [0, 2^31-1].

Mit diesem Programm habe ich den Code erzeugt:
static void Main(string[] args)
        {
            double U = 0xffffffff;
            double L = 0;
            double root;
            double p = 0.10;

            while (U > 1)
            {
                root = Math.Sqrt(U) * (1 - p);
                L = U * ((1 - p) * (1 - p)) / ((1 + p) * (1 + p));
                //Console.WriteLine("{0} - {1} = {2}", L, U - 1, root);
                Console.WriteLine("if (x >= {0}) return {2};", (uint)L, (uint)U - 1, (uint)root);
                U = L;
            }
           
            Console.ReadKey();
        }

Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: Svenska am 19. February 2017, 00:55
naja das iterationsverfahren muss schon 3 bis 4 mal durchgeführt werden
Naja, das Quake3-Verfahren (berechnet allerdings die Inverse und als Float) nutzt einen geschickten Anfangswert und kommt daher mit nur einer Iteration aus.

Das sind einige Maschinenbefehle die da durchgejagt werden müssen. Blöderweise auch Divisionen die relativ langsam sind. Das wird dann langsamer als gleich die FPU zu bemühen.
Auf modernen Maschinen ist die Anzahl der Instruktionen nicht unbedingt ein guter Indikator für die Gesamtgeschwindigkeit.

Wenn du den Definitionsbereich besser eingrenzen kannst, wäre vielleicht eine Tabelle (mit linerarer Interpolation) was für dich. Auf kleinen Mikrocontrollern erschlägt man die meisten komplexen Funktionen (sin, log, e^x) mit solchen Ansätzen.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 20. February 2017, 10:48
ja das Verfahren von John Carmack...

das mit einer Lut und einer linearen Interpolation funktioniert ganz gut bei Sin/Cos. da diese Funktionen relative linear sind.
Exponential /Wurzelfunktionen sind aber stark unter/über-Linear, und wenn man einen Eingangswertebereich von 32-Bit hat was mehr als 4 MRD. sind, wird es schwer. Nehmen wir mal an, wir nehmen nur jeden tausendsten Wert (immer noch 4 Mio-Werte) in die LUT auf
und interpolieren dazwischen, bekommt man rein gar nix mehr Vernünftiges. Man wird z.B. zwischen 0 und 1000 interpolieren.
Und eine LUT in der Größe ist auch nicht zu empfehlen weil sie das gesamte Caching durcheinanderbringt und in allen Ecken und Enden Cache-Misses erzeugt.

Ich brauche etwas, dass mir im Schnitt die Wurzel in 10 Takten berechnet, d.h. auf einem 3GHz Rechner 300 Mio. mal pro Sekunde.
mit den Voraussetzungen, 32Bit. Eingangsbandbreite mit (logischerweise) 16 Bit Ergebnisbandbreite.
und einer maximalen Abweichung von 10-15%. Wichtig ist vor allem, dass die Ergebnisse stetig sind, d.h. nicht irgendwo ein kleinerer Eingangswert ein größeres Ergebnis liefert als ein größerer Eingangswert.


 


Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: kevin am 20. February 2017, 11:59
Zehn Takte sind ja quasi gar nichts. Auf die nächste Zweierpotenz kriegst du damit wahrscheinlich hin (was eine Abweichung bis ungefähr 22,5% wäre, wenn ich das gerade richtig gerechnet habe), aber sonst sehe ich eher schwarz...
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 20. February 2017, 13:10
und wie sähe dann deine Rechnung aus mit etwa 22% Abweichung?
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: kevin am 20. February 2017, 17:22
Vergiss es, Denkfehler meinerseits. Du kriegst die Wurzel nicht von Zweierpotenzen, sondern nur von Quadraten von Zweierpotenzen, womit der Fehler nochmal ein gutes Stückchen größer wird.

Der Ansatz war, dass du mit bsr die Bitnummer des höchstwertigen gesetzten Bits kriegst. Bei den Quadraten von Zweierpotenzen ist die Wurzel dann (1 << (bitnr >> 1)), bei allem anderen geht es entsprechend mehr oder weniger daneben. Für "ungerade" Zahlen kriegst du n für alles bis (n + 1)² - 1, was natürlich ungenauer als nötig ist. Man will den Schritt irgendwo in der Mitte zwischen zwei Quadraten von Zweierpotenzen machen, was durch passende Rundung machbar sein sollte.

Mit nur zehn Takten wäre allerdings schon das eher eng, glaube ich...
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: xenos1984 am 20. February 2017, 18:04
Wichtig ist vor allem, dass die Ergebnisse stetig sind, d.h. nicht irgendwo ein kleinerer Eingangswert ein größeres Ergebnis liefert als ein größerer Eingangswert.
Diese Eigenschaft nennt man nicht stetig, sondern monoton ;) Aber das nur als Randbemerkung.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: Svenska am 21. February 2017, 02:46
Exponential /Wurzelfunktionen sind aber stark unter/über-Linear, und wenn man einen Eingangswertebereich von 32-Bit hat was mehr als 4 MRD. sind, wird es schwer
Darum schrieb ich, dass du nach Möglichkeit den Definitionsbereich besser eingrenzen solltest.

Nehmen wir mal an, wir nehmen nur jeden tausendsten Wert (immer noch 4 Mio-Werte) in die LUT auf und interpolieren dazwischen, bekommt man rein gar nix mehr Vernünftiges.
Du musst die Stützpunkte ja nicht gleichmäßig über den gesamten Definitionsbereich verteilen. Oder du nimmst zwei Tabellen fixer Größe, eine für Zahlen bis 10000 und eine für Zahlen darüber (wo die Nichtlinearität nicht mehr so weh tut).

Es hängt vom Anwendungsfall ab, was praktikabel ist und was nicht, aber darüber schweigst du dich ja netterweise aus...

Ich brauche etwas, dass mir im Schnitt die Wurzel in 10 Takten berechnet, d.h. auf einem 3GHz Rechner 300 Mio. mal pro Sekunde.
Das ist reichlich sportlich. Ich rate dir einfach mal, dein Design grundsätzlich darauf abzuklopfen, ob du das Problem nicht anderweitig lösen kannst.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 21. February 2017, 11:09
Für was das gut sein soll?
hat ja keiner danach gefragt. :roll:
Also es geht um das erkennen von komplexen Datenkorrelationen in einem KI-System. Und speziell hier um das Erkennen von unter/über-linearen Abhängigkeiten von Daten. (Das mal vorweg) :-)

Die FPU braucht ca. 50-60 Takte für eine Wurzelberechnung. Dafür ist das Ergebnis exakt und auch richtig gerundet.
und bei mir geht's nicht um Genauigkeit, da es eher um das Ermitteln von Wahrscheinlichkeiten geht.

Die geteilte LUT - Funktion wäre die schneller als 50 Takte? Sie müsste ja auch bedingte Sprünge enthalten und Zugriffe auf den Arbeitsspeicher. Die Interpolation benötigt auch ein paar Rechenoperation auch mul und div, was nicht gerade schnell ist.
Oder täusche ich mich da?
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: MasterLee am 21. February 2017, 11:19
Ich komme leider nicht unter 81 Befehle ;(
Dafür ist es aber exakt
Eingabe R0 Ausgabe R1. R2 Wird während der Berechnung überschrieben

;Zu berechnendes X
MOV R0,#4096
;bit= #0x40000000
MOV R1,#0 ;res

ADD R2,R1,#0x40000000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x40000000

ADD R2,R1,#0x10000000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x10000000

ADD R2,R1,#0x4000000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x4000000

ADD R2,R1,#0x1000000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x1000000

ADD R2,R1,#0x400000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x400000

ADD R2,R1,#0x100000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x100000

ADD R2,R1,#0x40000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x40000

ADD R2,R1,#0x10000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x10000

ADD R2,R1,#0x4000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x4000

ADD R2,R1,#0x1000
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x1000

ADD R2,R1,#0x400
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x400

ADD R2,R1,#0x100
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x100

ADD R2,R1,#0x40
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x40

ADD R2,R1,#0x10
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x10

ADD R2,R1,#0x4
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x4

ADD R2,R1,#0x1
SUBS R2,R0,R2
MOV R1,R1,LSR #1
MOVpl R0,R2
ADDpl R1,R1,#0x1

Aber wenn man das ganze in VHDL rekodiert kann man es bestimmt in einen Takt schaffen…
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: Svenska am 23. February 2017, 04:49
Also es geht um das erkennen von komplexen Datenkorrelationen in einem KI-System.
Bleibt die Frage, ob du den Definitionsbereich irgendwie eingrenzen kannst, bzw. wo die Genauigkeitsanforderungen am höchsten sind. Deine Erklärung hilft mir da nicht weiter. :-D

Wenn du z.B. weißt, dass deine Werte fast immer klein sind, dann darf die Funktion für große Zahlen auch extreme Werte liefern. Oder du nimmst für große Zahlen die FPU und nur für die kleinen Zahlen eine schnelle Abschätzung per LUT. Ob das für dich praktikabel ist, kann ich nicht einschätzen.

Die geteilte LUT - Funktion wäre die schneller als 50 Takte? Sie müsste ja auch bedingte Sprünge enthalten und Zugriffe auf den Arbeitsspeicher.
Keine Ahnung, ob das schneller ist. Mit x86 habe ich nicht so viel zu tun. Mit Sprungvorhersage und einer kleinen, dauerhaft im Cache liegenden LUT sollte das aber nicht weh tun.

Die Interpolation benötigt auch ein paar Rechenoperation auch mul und div, was nicht gerade schnell ist.
Nö, das sind nur wenige Befehle (in Pseudo-Assembler, ungetestet):
; Ein-/Ausgabe in R1; clobbers R2, R3, R4
LOOKUP:
  MOV R4, #LUT      ; Pointer zur LUT
  SHR R1, R1, #24   ; obere 8 Bit der Eingabe als Tabellenindex (abgerundet)
  LD R2, [R4 + R1]  ; lade a
  INC R1            ; nächster Index
  LD R3, [R4 + R1]  ; lade b
  ADD R1, R2, R3    ; linerare Interpolation ist "(a+b) / 2"
  SHR R1, #1        ; kein DIV nötig

Wenn dich die Nichtlinearität für kleine Zahlen stört, musst du für die Indizierung etwas mehr Gehirnschmalz investieren, indem du nicht gleichmäßig verteilst. Da könnte man z.B. mit BSR (Bit Scan Reverse == log base 2) oder CLZ (Count Leading Zeros) arbeiten, um den Index auszurechnen. Bei sowas ist z.B. Excel dein Freund, damit kannst du sowohl die Genauigkeit abschätzen als auch die Tabelle selbst mit generieren.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 23. February 2017, 11:05
Lineare Interpolation ist nicht: (a+b)/2   !!

Lineare Interpolation berücksichtigt wie nahe der Wert an den jeweiligen Grenzen liegt.
Und je nachdem wird er dann gewichtet.
Darum ist da eine Division mit drin und paar Subraktionen.

Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: MasterLee am 23. February 2017, 18:15
Ist lineare interpolation nicht:
(1-t)*a+t*bbzw.
a+t*(b-a)Wo siehst du da eine division ich sehe nur eine Multiplikationen?
Bzw. wenn man mit Ganzahlen arbeitet und t von 0 bis 256 geht:
((a<<8)+t*(b-a))>>8
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: Svenska am 23. February 2017, 18:25
Stimmt, ich dachte an den Mittelwert. War halt spät. :roll:
Die relevante Frage ist, ob das die Anforderungen an Genauigkeit und Geschwindigkeit erfüllt, und das müsste man austesten.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: MasterLee am 23. February 2017, 19:51
Stimmt, ich dachte an den Mittelwert. War halt spät. :roll:
Die relevante Frage ist, ob das die Anforderungen an Genauigkeit und Geschwindigkeit erfüllt, und das müsste man austesten.
ganz falsch ist dein Ansatz ja nicht mit wiederholter Ermittelung des Mittelwerts lässt sich auch Interpolieren man fängt dann beim höchsten Bit an und arbeitet sich nach unten durch ist ein Bit gesetzt nimmt man:
a'=Mittewert(a+b)
b'=b
wenn es nicht gesetzt ist nimmt man
a'=a
b'=Mittewert(a+b)
Am ende enthält a das Ergebnis, das das selbe wäre wie bei:
((a<<8)+t*(b-a))>>8Wenn t 8 Bit groß ist und obige Methode somit 8 mal durchgeführt würde.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 24. February 2017, 13:31
hmmm, ja gut aber alles nicht wirklich schneller also die 50 Takte der FPU.

10 Takte auf einem IA-32 (x86) sollten doch irgendwie möglich sein...

Oder gehört das ins Reich der Träume. :-D
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: OsDevNewbie am 25. February 2017, 23:12
Du könntest die x86-Instruktion "sqrtss" verwenden. Sie berechnet eine Quadratwurzel aus einer 32-bit Floating-Point-Zahl und braucht so 13 Zyklen auf modernen Prozessoren. Mit "sqrtps" kannst du sogar 4 Zahlen gleichzeitig berechnen in der gleichen Zeit. Das führt dazu, dass du 13/4 = 3.25, also etwa 4 Zyklen pro Zahl brauchst. Natürlich geht das nur, wenn du SSE verwenden kannst.
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: Svenska am 26. February 2017, 17:20
Hast du gemessen oder geschätzt?
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: TomCat am 27. February 2017, 09:11
Hast du da auch den Code dazu?
Wäre sehr interessant. :-)
Titel: Re: Quadratwurzel schnell berechnen
Beitrag von: OsDevNewbie am 01. March 2017, 14:12
Hast du gemessen oder geschätzt?
Ich habe im Intel optimization manual nachgeschaut. Dort sind aber leider nur die Latenzen auf Prozessoren der letzten paar Generationen verfügbar.

Hast du da auch den Code dazu?
Wäre sehr interessant. :-)
Ich habe im Anhang ein C Programm (unkommentiert, aber hoffentlich leicht lesbar), dass auf verschiedene Arten die Wurzel aus einer int Zahl zieht. Bei meiner Testmaschine (i7-4770k mit 32GiB RAM) kommen folgende Ergebnisse raus:
Alloziere 3000000000 (11718750.00 KiB) Elemente...OK
Initialisiere...OK
Berechne (calc_sqrt_float_sse)...OK(1.642945 Zyklen/Element)
OK
Berechne (calc_sqrt_float_avx)...OK(1.619843 Zyklen/Element)
OK
Berechne (calc_sqrt_double_sse)...OK(6.353434 Zyklen/Element)
OK
Berechne (calc_sqrt_double_avx)...OK(6.327483 Zyklen/Element)
OK
Berechne (calc_sqrt_float_sse_2x)...OK(1.621323 Zyklen/Element)
OK
Berechne (calc_sqrt_float_avx_2x)...OK(1.625151 Zyklen/Element)
OK
Berechne (calc_sqrt_double_sse_2x)...OK(6.315364 Zyklen/Element)
OK
Berechne (calc_sqrt_double_avx_2x)...OK(6.325675 Zyklen/Element)
Wie du siehst ist es sehr schnell. (Damit die AVX-Funktionen ausgeführt werden musst du beim Compilieren AVX aktivieren).
(Wieso die AVX-Variante gleich schnell ist wie die SSE-Variante verstehe ich selber noch nicht ganz.)
Wenn man die Wurzel aus floats zieht, dann hat man eine Genauigkeit von 24-Bit ich weiss nicht ob dir das reicht. Deshalb hab ich noch die double Varianten davon gemacht, da hat man eine höhere Genauigkeit (wie hoch steht nicht im Manual), aber sind teilweise ein bisschen langsamer.
Du kannst das Programm ja mal bei dir testen und schauen, ob es schnell genug ist für dich.
Du darfst auch von meinem Code kopieren oder mir sagen, wenn ich etwas gemacht habe.

Ich muss noch dazu sagen, dass ich jetzt die Umwandlung zwischen ints und floats in meinen Berechnungen in meinem Post nicht einberechnet habe (die kosten immer so 4 Zyklen pro umwandlung), sie werden aber im Programm mit eingerechnet.