Quadratwurzelberechnung
Rainer Glaschick
2011-10-15, 2022-09-02
Die Berechnung der Quadratwurzel ist bekannt; dieser Text fasst lediglich meine Sicht zusammen und behandelt die Umsetzung auf Digitalrechnern.
Mit `a` oder `b^2` wird die Zahl bezeichnet, deren Quadratwurzel bestimmt werden soll.
Inhalt:
1. Iteratives Verfahren nach Heron
Dies Verfahren ist der im folgenden Kapitel beschriebenen Nullstellenbestimmung nach Newton äquivalent.
Heron-Folge
Die Folge
- `x_(i+1) = 1/2 * (x_i + b^2 / x_i )`
hat den Fixpunkt `x=b=sqrt(a)` und konvergiert gegen den Fixpunkt.
Im folgenden wird `y` anstelle von `x_(i+1)` geschrieben und dann der Index bei `x` weggelassen:
- `y = 1/2 * (x + b^2 / x) = (x^2 + b^2)/(2*x)`
Zum Beweis wird die Differenz `z = x - y` zweier Werte betrachtet:
- `z = x - x/2 - b^2 / (2*x) = (x - b^2/x) / 2`
- `z = (x^2 - b^2) / (2*x) = ((x+b)*(x-b))/(2*x)`
Ersichtlich verschwindet die Differenz bei `x=b` (Fixpunkt).
Sofern `x>b` ist, ist auch `z>0`; d.h., wenn der aktuelle Schätzwert größer als der gesuchte Wert ist, dann wird er kleiner und bleibt nicht gleich; der Folgewert ist also monoton fallend.
Die Folge ist für `x>b` auch durch `b` beschränkt:
- `y - b = (x^2 + b^2 -2*b*x) / (2*x) = ((x-b)^2)/(2*x) > 0`
Da die Differenz immer (außer bei `x=b`) positiv ist, liegen alle Folgenglieder außer der ersten immer oberhalb der Wurzel; die Folge konvertiert also immer.
Im folgenden wird daher meist impliziert, dass `x>b` ist.
Sofern der Wert von `x` mehr als doppelt so groß wie der Fixpunkt liegt, wird:
- `y ~~ x^2 / (2x) = x/2`
D.h., der Wert wird mit jedem Schritt halbiert, bis er in die Nähe des Fixpunktes kommt. In dessen Nähe konvergiert die Folge dann erheblich schneller.
Für den Startwert `x=1` und wesentliche größeres Argument `b^2` ist der nächste Wert:
- `y = (1 + b^2)/2 ~~ b^2 / 2` ist.
Das entspricht dem Startwert `x=b^2`. Es ist also im allgemeinen Fall für die Konvergenz gleich, ob mit dem Startwert `x=1`, `x=a` oder `x=½*a` begonnen wird. Mit `x=1` wird ein geringer Vorteil für kleine Argumente bewirkt; mit `x=a` einer für große Argumente. Der Startwert `x=½*a` beginnt dann mit `1<a<2` mit einem Anfangswert kleiner als eins; dieser wird aber im ersten Durchlauf ja sofort größer als Eins.
Allerdings hat der Anfangswert x=a den Vorteil, dass — wenn a>1 — die Werte von x monoton fallen, solange die Rechengenauigkeit ausreicht;
fallen sie nicht mehr monoton, kann die Iteration abgebrochen werden.
Die Funktion sieht dann im Kern so aus,
wobei eps die gewünschte (absolute) Genauigkeit ist:
sqrt(a):
if a < 1
return (1 / sqrt(1/a))
x = a
while x - a/x > eps
x = (x + a/x) / 2
return x
Die meisten modernen Compiler optimieren die doppelte Berechnung
von a/x; ansonsten würde man schreiben:
sqrt(a):
if a < 1
return 1/sqrt(1/a)
x = a
t = 1
while x - t > eps
t = a/x
x = (x + t) /2
return x
Hier sieht man auch, dass nicht eps=0 gesetzt werden kann,
da sich x und t durch lediglich
ein Bit in der Mantisse unterscheiden können,
welches bei der Division durch zwei wieder verloren geht,
so dass x immer größer als t bleibt.
Hingegen kann man als Abbruchkriterium verwenden,
dass die Werte für x monoton fallen;
tauchen zwei gleiche x auf, ist anzuhalten:
sqrt(a):
if a < 1
return 1/sqrt(1/a)
x = a
xx = x+1
while x < xx
xx = x
x = (x + a/x) /2
return x
Hier werden in der Schleife vier Fließkommaberechnungen (Vergleich eingeschlossen) pro Umlauf benötigt.
Da bei Fließkommaarithmetik die Multiplikation nicht viel aufwendiger ist als eine Addition, kann anstelle eines additiven Deltas auch mit einem Verbesserungsfaktor gearbeitet werden:
- `y = v*x`
- `v = y/x = (x^2 + b^2)/(2*x^2) = (1 + b^2/x^2)/2`
Wie vor gilt, dass, da der Faktor für `b=x` Eins wird, die Wurzel ein Fixpunkt ist. Weiterhin ist für `x>b` der Faktor `v<1`.
Die Funktion sieht dann so aus:
sqrt(a):
v = 0; x = a
while v < 1
v = (1 + a/(x*x)) / 2
x *= v
return x
Die Abfrage v<1 müsste eher v < 1-epsilon lauten;
man kann aber alternativ auch den vorigen Wert mithalten
und auf monotone Vergrößerung abfragen:
sqrt(a):
v = 0
vv = 1
x = a
while v < vv
vv = v
v = (1 + a/(x*x)) / 2
x *= v
return x
Besser (für gerundete Fließkommaarithmetik):
sqrt(x):
x = a
xx = a+1
while xx > x
xx = x
x = (x + a/x) /2
return x
Mit Genauigkeitsangabe eps (bezogen auf 1):
sqrt(x, eps):
eps1 = 1 + eps
x = a
xx = x * eps1 + eps
while xx > x*eps1
xx = x
x = (x + a/x) / 2
return x
Anfangswertbestimmung
Ist der Zahlenbereich bekannt, dann ist auch der größte sinnvolle Anfangswert gegeben, z.B. bei `2^32` ist es `2^16`.
Ein mögliche Lösung wäre dann eine Tabelle der Potenzen von zwei und ihrer Wurzeln.
Statt dessen könnten die Werte aufgezählt werden:
sqrt(a):
x = 1
while a/x > x: x *= 2
Hier nur drei anstelle von vier Operationen benötigt, aber die Anzahl der Durchläufe liegt in der gleichen Größenordnung. Gleich, ob der Näherungswert der Wurzel durch Verdoppeln oder Halbieren erreicht wird, es werden immer ungefähr `½*ld(a)` Schritte benötigt.
Vorskalierung
Anstelle einer Bestimmung des Anfangswerts ist es sehr effizient, zunächst mit einer als Code realisierten kleinen Tabelle die Berechnung auf kleinere Zahlen zu reduzieren:
sqrt(a):
if a < 1 : return 1/(sqrt 1/a)
if a < 100 : return 10*(sqrt a/100)
if a < 10000 : return 100*(sqrt a/10000)
if a < 1000000: return 1000*(sqrt a/10000000
\ weiter wie oben
Damit wird z.B. bei Zahlen über 1 Million die Anzahl der Durchläufe von ca. 26 auf 6 reduziert. Die Tabelle braucht nicht weiter ausgebaut zu werden, da die Genauigkeit der Fließkommaarithmetik ohnehin die Anzahl der Durchläufe begrenzt: Bei einer Mantisse von 52 bit (double) ist die Anzahl auf 53 beschränkt. Da `10^6` ungefähr `2^20` ist, sind weitere Verbesserungen nicht mehr möglich, wenn die Mantisse 23 Bit lang ist (single); am besten wären hier wohl zwei Abfragen mit `2^8` und `2^16`. Bei einer 52 Bit Mantisse sind das `2^12`, `2^24` und `2^36`.
Rechengenauigkeit und Überlauf
2. Nullstellenbestimmung nach Newton
Das Verfahren nach Newton (und Raphson) bestimmt iterativ die Nullstelle einer differenzierbaren Funktion `f(x)`, indem der aktuelle Schätzwert durch den Quotienten aus dem aktuellen Funktionswert und der Ableitung an diesem Punkt vermindert wird:
- `x_(n+1) = x_n - f(x) / (f'(x))
Erste Formel
Gesucht wird die Nullstelle von
- `f(x) = x^2 - a`
Die Ableitung ist
- `f'(x) = 2x
Damit ist die abzuziehende Verbesserung:
- `y_n = (f(x_n) / (f'(x_n))) = (x_n^2 -a) / (2x_n)
und der neue Wert
- `x_(n+1) = x_n - (x_n^2 - a) / (2*x_n) = (x_n^2 + a) / (2x_n) = 1/2 (x_n + a/x_n)
und ist damit dem iterativen Verfahren von Heron gleich.
Zweite Formel
Gesucht wird die Nullstelle von
- `f(x) = x - a/x
Die Ableitung ist
- `f'(x) = 1 + a/(x^2)
Sinnvollerweise ist `x>=1` und `a >= 1`. Der verbesserte Wert `y` ist:
- `y = x - (f(x) / (f'(x))) = x - x (x^2 - a) / (x^2 + a) = x (1 - (x^2-a) / (x^2+a))
Dritte Formel
Gesucht wird die Nullstelle von
- `f(x) = a/x^2 - 1`
Die Ableitung ist
- `f'(x) = -2*a/x^3`
Damit ist der verbesserte Wert
- `y = x - (f(x) / f'(x)) = x - (a/x^2 -1) / (-2*a/x^3)`
- `y = x + (a*x - x^3)/(2*a) = x + x/2 - x^3/(2*a)`
- `y = x + x/2 (1 - 1/a x^2 )
Wenn man den Kehrwert von `a` vorausberechnet, dann werden pro Schritt eine Quadrierung und zwei Multiplikationen, aber keine Divisionen benötigt.
Das Delta `z` ist
- `z = y - x = x/2 (1 - x^2 / a)
- `x^2/a = 1 => z = 0: "Fixpunkt"
- `x^2/a > 1 <=> x^2 > a => z < 0: "fallend"
- `x^2/a < 1 <=> x^2 < a => z > 0: "steigend"
Für den Fall `a > 1` wird durch `x = 1` die Bedingung `x^2 < a` erfüllt und eine monoton steigende Folge erreicht, solange `x^2 < a` ist, d.h. solange die Wurzel nicht gefunden ist; Abbruchkriterium ist entweder `z < epsilon` oder besser `x^2 >= a`, weil so kein `epsilon` benötigt wird. Eine Abfrage auf Gleichheit ist problematisch, da durch Rundungsfehler nicht auszuschließen ist, dass `x^2` geringfügig größer ist als `a`.
Im anderen Fall `a < 1` ist dies nicht so einfach, weil `x=0` entfällt; aber hier kann die Wurzel vom Kehrwert bestimmt und davon der Kehrwert genommen werden.
Es ist `y` durch `b` nach oben beschränkt:
- `b - y = b - x - x/2 ( 1 - x^2/b^2) = (b-x) -x/2(b^2 -x^2)/(b^2) = (b-x) (1-(bx +x^2)/(2b^2))
- `x < b => (bx+x^2)/(2b^2) < 1 => b-y > 0
Wegen `b = sqrt(a)` ist die Folge nach oben beschränkt, solange `x < sqrt(a)` ist; d.h. solange der Anfangswert diese Bedingung erfüllt, konvergiert die Folge gegen `b = sqrt(a)`.
Zur Bestimmung der Konvergenz:
- `(b-y)/(b-x) = (2b^2-bx-x^2)/(2b^2)
Sei `beta = b-x; x = b - beta`, dann ist
- `(b-y)/(b-x)^2 = (2b^2-b(b-beta) - (b-beta)^2)/(2beta b^2) = (3 beta b - beta^2) / ( 2 beta b^2) = (3b -beta) /(2b^2)
- `(b-y)/(b-x)^2 ~~ 3/2 cdot 1/b if beta << b
Also konvergiert die Folge quadratisch.
3. Berechnung durch sukzessive Approximation
Diese Methode wurde bereits von Konrad Zuse bei dem Rechner Z3 verwendet; sie kann nützlich sein, wenn eine Division nicht vorhanden oder aufwändig ist.
Es werden die Potenzen von 2 absteigend verprobt, ob ihre Addition das Ergebnis verbessert:
sqrt(n):
assert n<2^31
b = 2^15
r = 0
while b > 0
if (r+b)^2 <= n
r =+ b
b = b / 2
return r
Ein aufsteigendes Verproben ist nicht möglich, wie man sich am Beispiel der Wurzel aus 16 leicht klarmacht.
Hingegen kann vorab die größte Potenz von 2 ermittelt werden, deren Quadrat kleiner ist:
int sqrt(n):
b = 2
while b^2 < n
b = b * 2
r = 0
while b > 0
if (r+b)^2 <= n
r =+ b
b = b / 2
return r
Dies kann aber zum Überlauf bei der Berechnung von b^2 führen,
wenn n nahe am Ende des Rechenbereichs liegt;
dies kann vermieden werden, ist aber geringfügig langsamer:
int sqrt(n):
b = 2
while b^2 <= n
if b^2 = n
:> b
b = b * 2
r = 0
while b > 0
if (r+b)^2 <= n
r =+ b
b = b / 2
return r
Dies kann leicht für Wurzeln höheren Grades erweiteret werden:
int (e) root (n):
b = 2
while b^e <= n
b = b * 2
r = 0
while b > 0
if (r+b)^e <= n
r =+ b
b = b / 2
return r
Bei einer Fließkomma-Arithmetik mit Zweier-Exponenten ist gleichfalls eine Multiplikation mit 2 bzw. Division durch 2 sehr einfach.
Da die Heron-Folge für die Quadratwurzel weniger Iterationen benötigt als die Intervallschachtelung, ist sie effzienter, wenn eine Division verfügbar ist.
4. Iterative Akkumulation
Für Ganzzahlen und Festkommaarithmetik kann die Quadratwurzel ohne Multiplizieren berechnet werden, da die Summe der ungeraden Zahlen das Quadrat der Anzahl ist:
- `sum_(k=1)^n (2k-1) = n^2`
- `"weil " (n+1)^2 = n^2 + 2n + 1 = n^2 + (2(n+1) -1)`
Anstatt die ungeraden Zahlen aufzusummieren
(und um Probleme mit dem Rechenbereich zu vermeiden),
werden sie von dem Argument x abgezogen, solange der Rest groß genug ist:
y = 1
z = 1
while x > y
x = x - y
y = y + 2
z = z + 1
Hier wird, wenn 'x' keine Quadratzahl ist, die nächst größere Zahl bestimmt.
Man kann das Register für die nächste ungerade Zahl sparen,
da die nächste ungerade Zahl 2*n-1 ist,
und daher zweimal n abgezogen und 1 addiert werden kann:
y = 1
while x > 0
x = x - y - y + 1
y = y + 1
Da die Berechnung `O(sqrt n)` Additionen benötigt, ist sie im Normalfall zu ineffizient.
Auf Basis dieser Methode, jedoch erheblich optimiert, wurde in der ENIAC die Wurzelberechnung realisiert: https://bshelburne.wittenberguniversity.org/HowTheENIACTookASquareRoot011909.pdf
5. Wurzel höheren Grades
Zur Berechnung der n-ten Wurzel wird wieder die Newton-Iteration verwendet:
- `f(x) = x^n - a
- `f'(x) = n x^(n-1)
- `x_(n+1) = x_n - (x^n -a) / (n x^(n-1)) = x_n - 1/n (x_n - a/x^(n-1)) = 1/n ((n-1)x_n + a / ( x_n^(n-1)))
