Wie Hypot2 (x, y) Berechnung, wenn Zahlen überlaufen können

Question

Ich möchte eine Hypot2-Berechnung auf einem 16-Bit-Prozessor tun.

Die Standardformel lautet c = sqrt((a * a) + (b * b)). Das Problem dabei ist, dass es bei großen Eingaben überläuft. Z.B. 200 und 250, multipliziere 200 * 200, um 90.000 zu erhalten, was höher ist als der maximale signierte Wert von 32.767, so dass es überläuft, wie auch b, werden die Zahlen addiert und das Ergebnis kann ebenso nutzlos sein; es könnte sogar eine Fehlerbedingung wegen eines negativen sqrt signalisieren.

In meinem Fall habe ich mit 32-Bit-Zahlen zu tun, aber 32-Bit-Multiplikation auf meinem Prozessor ist sehr schnell, etwa 4 Zyklen. Ich benutze einen dsPIC Mikrocontroller. Ich möchte lieber nicht mit 64-Bit-Zahlen multiplizieren, weil das wertvolles Gedächtnis verschwendet und zweifellos langsamer sein wird. Außerdem habe ich nur sqrt für 32-Bit-Nummern, so dass 64-Bit-Nummern eine andere Funktion benötigen. Also, wie kann ich einen Hypot berechnen, wenn die Werte groß sein können?

Bitte beachten Sie, dass ich nur ganzzahlige Mathematik dafür verwenden kann. Die Verwendung von Fließkomma-Mathe führt zu einem Geschwindigkeitshit, den ich eher vermeiden möchte. Mein Prozessor hat eine schnelle Integer/Fixed-Point atan2-Routine, etwa 130 Zyklen; Könnte ich damit die Länge der Hypotenuse berechnen?

Answer 1

Je nachdem, wie viel Genauigkeit Sie benötigen, können Sie die Quadrate und die Quadratwurzeloperation vermeiden. Es gibt einen Abschnitt zu diesem Thema in (Abschnitt 10.2, "High-Speed ​​Vector-Magnitude Approximation", ab Seite 400 in meiner Ausgabe).

Die Näherung ist im wesentlichen:

magnitude = alpha * min + beta * max 

wobei max und min die Maximal- und Minimalabsolutwerte der realen und imaginären Komponenten sind, und alpha und beta sind zwei Konstanten, die ausgewählt werden, eine angemessene Fehlerverteilung zu geben, über den Bereich von Interesse.Diese Konstanten können als Brüche mit einer Potenz von 2 Teilern dargestellt werden, um die Arithmetik einfach/effizient zu halten. In dem Buch schlägt er vor alpha = 15/16, beta = 15/32, und Sie können dann die Formel vereinfachen:

magnitude = (15/16) * (max + min/2) 

die umgesetzt werden könnten wie folgt Integer-Operationen verwenden:

magnitude = 15 * (max + min/2)/16 

und natürlich können wir Schichten für die Divisionen verwenden:

magnitude = (15 * (max + (min >> 1))) >> 4 

Fehler ist +/- 5% über einem Quadranten.

Weitere Informationen zu dieser Technik hier: http://www.dspguru.com/dsp/tricks/magnitude-estimator

Answer 2

Benötigen Sie volle Präzision? Wenn Sie dies nicht tun, können Sie Ihre Reichweite ein wenig erhöhen, indem Sie einige weniger signifikante Bits verwerfen und danach multiplizieren.

Können a und b alles sein? Wie wäre es mit einer Nachschlagetabelle, wenn Sie nur ein paar a und b haben, die Sie berechnen müssen?

Answer 3

Eine einfache Lösung Überlauf zu vermeiden, ist sowohl a und b durch a+b vor Quadrierung zu teilen und dann durch die Quadratwurzel a+b multiplizieren. Oder machen Sie das gleiche mit max(a,b).

Answer 4

Sie können ein wenig einfache Algebra tun, um die Ergebnisse wieder in den Bereich zu bringen.

sqrt((a * a) + (b * b)) 
= 2 * sqrt(((a * a) + (b * b))/4) 
= 2 * sqrt((a * a)/4 + (b * b)/4) 
= 2 * sqrt((a/2 * a/2) + (b/2 * b/2)) 

Answer 5

Aniko und John, es scheint mir, dass Sie das OP-Problem nicht angesprochen haben. Wenn a und b ganze Zahlen sind, ist a * a + b * b wahrscheinlich überlaufen, weil Integer-Operationen ausgeführt werden. Die offensichtliche Lösung besteht darin, a und b in Gleitkommawerte zu konvertieren, bevor a * a + b * b berechnet wird. Aber das OP hat uns nicht wissen lassen, welche Sprache wir verwenden sollten, also sind wir ein bisschen festgefahren.

Answer 6

Da Sie im Grunde keine Multiplikationen ohne Überlauf tun können, werden Sie wahrscheinlich etwas Genauigkeit verlieren.

Um die Zahlen in einem akzeptablen Bereich zu bekommen, x einige Faktor herausziehen und

c = x*sqrt((a/x)*(a/x) + (b/x)*(b/x)) 

verwenden Wenn x ein gemeinsamer Faktor ist, werden Sie nicht Präzision verlieren, aber wenn es nicht, dann ist werden Sie Präzision verlieren .

Update: Noch besser ist, da Sie einige leichte Arbeit mit 64-Bit-Zahlen zu tun, mit nur einem 64-Bit Außerdem kann man ein winzigen mit nur den Rest dieses Problems in 32-Bit zu tun Verlust der Genauigkeit. Um dies zu tun: Machen Sie die zwei 32-Bit-Multiplikationen, um Ihnen zwei 64-Bit-Zahlen zu geben, fügen Sie diese und dann Bit-Verschiebung wie benötigt, um die Summe zurück auf 32 Bits vor der Wurzel zu nehmen. Wenn Sie Bit immer um 2 Bits verschieben, multiplizieren Sie einfach das Endergebnis mit 2^(die Hälfte der Bitverschiebungen), basierend auf der obigen Regel.Das Abschneiden sollte nur einen sehr geringen Genauigkeitsverlust verursachen, nicht mehr als 2^31 oder 0,00000005% Fehler.

Answer 7

Dies ist wörtlich aus this @John D. Cook blog post, daher CW genommen:

Hier ist, wie sqrt(x*x + y*y) ohne Überlauf zu riskieren, zu berechnen.

1. max = maximum(|x|, |y|)
2. min = minimum(|x|, |y|)
3. r = min/max
4. return max*sqrt(1 + r*r)

Wenn @John D. Koch entlang und Pfosten kommt dies sollte man ihm geben das akzeptieren :)

Answer 8

Die Standardformel ist c = sqrt ((a * a) + (b * b)). Das Problem liegt bei großen> Eingängen, die überlaufen.

Die Lösung für Überläufe (abgesehen von einem Fehler) ist, Ihre Zwischenberechnungen zu sättigen.

Berechnen Sie C = a * a + b * b. Wenn a und b 16-Bit-Nummern sind, haben Sie nie einen Überlauf. Wenn es sich um vorzeichenlose Nummern handelt, müssen Sie die Eingaben zuerst nach rechts verschieben, damit die Summe in eine 32-Bit-Zahl passt.

Wenn C> (MAX_RADIUS)^2, MAX_RADIUS zurückgeben, wobei MAX_RADIUS der Maximalwert ist, den Sie tolerieren können, bevor Sie einen Überlauf auslösen.

Ansonsten verwenden Sie entweder sqrt() oder CORDIC algorithm, was die Kosten von Quadratwurzeln zugunsten von Schleifeniteration vermeidet und + Verschiebungen addiert, um die Amplitude des (a, b) Vektors abzurufen.

Answer 9

Wenn Sie a und b auf höchstens 7 Bits beschränken können, werden Sie keinen Überlauf erhalten. Sie können eine Zählung-führende-Nullen-Anweisung verwenden, um herauszufinden, wie viele Bits weggeworfen werden sollen.

Angenommen a> = b.

int bits = 16 - count_leading_zeros(a); 
if (bits > 7) { 
    a >>= bits - 7; 
    b >>= bits - 7; 
} 
c = sqrt(a*a + b*b); 
if (bits > 7) { 
    c <<= bits - 7; 
} 

Viele Prozessoren haben diese Anweisung heutzutage, und wenn nicht, können Sie andere fast techniques verwenden.

Obwohl dies nicht die genaue Antwort geben wird, wird es sehr nahe (höchstens ~ 1% niedrig) sein.