Retro Computer Laboratory

Diese Seite ist nicht editierbar. Sie können den Quelltext sehen, jedoch nicht verändern. Kontaktieren Sie den Administrator, wenn Sie glauben, dass hier ein Fehler vorliegt.
<nimla>
Quadratwurzelberechnung

Autor: Rainer Glaschick
Datum: 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:

\Inhalt

1. Iteratives Verfahren nach Heron
======

Dies Verfahren ist der im folgenden Kapitel beschriebenen [Nullstellen...]
ä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=&frac12;*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=&frac12;*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 `&frac12;*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
-------


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








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.

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]
     

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



\ASCIIMATHML ASCIIMathML.js

</nimla>