Thema: Quadratwurzel schnell berechnen

[TomCat]

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

[kevin]

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.

[TomCat]

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.

 

[kevin]

Wenn du einen Float kriegen kannst (und dann auf der Binärdarstellung davon mit Integeroperationen arbeiten), ist eventuell das hier was für dich. Soll dann am Ende irgendwie so aussehen:

Code: [Auswählen]

val_int = (1 << 29) + (val_int >> 1) - (1 << 22) + a;

[Jidder]

Ich hätte da auch noch einen Vorschlag (für Fehler < 10%):

Code: [Auswählen]

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:

Code: [Auswählen]

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();
        }

[Svenska]

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.

[TomCat]

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.


 

[kevin]

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...

[TomCat]

und wie sähe dann deine Rechnung aus mit etwa 22% Abweichung?

[kevin]

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...

[xenos1984]

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.

[Svenska]

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.

[TomCat]

Für was das gut sein soll?
hat ja keiner danach gefragt.
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?

[MasterLee]

Ich komme leider nicht unter 81 Befehle ;(
Dafür ist es aber exakt
Eingabe R0 Ausgabe R1. R2 Wird während der Berechnung überschrieben

Code: [Auswählen]

;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…

[Svenska]

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.

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

Code: [Auswählen]

; 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.

[TomCat]

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.

[MasterLee]

Ist lineare interpolation nicht:

Code: [Auswählen]

(1-t)*a+t*bbzw.

Code: [Auswählen]

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:

Code: [Auswählen]

((a<<8)+t*(b-a))>>8

[Svenska]

Stimmt, ich dachte an den Mittelwert. War halt spät.
Die relevante Frage ist, ob das die Anforderungen an Genauigkeit und Geschwindigkeit erfüllt, und das müsste man austesten.

[MasterLee]

Stimmt, ich dachte an den Mittelwert. War halt spät.
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:

Code: [Auswählen]

a'=Mittewert(a+b)
b'=bwenn es nicht gesetzt ist nimmt man

Code: [Auswählen]

a'=a
b'=Mittewert(a+b)Am ende enthält a das Ergebnis, das das selbe wäre wie bei:

Code: [Auswählen]

((a<<8)+t*(b-a))>>8Wenn t 8 Bit groß ist und obige Methode somit 8 mal durchgeführt würde.

[TomCat]

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.