Heron-Verfahren

Das Babylonische Wurzelziehen (oft auch Heron-Verfahren) ist ein alter iterativer Algorithmus zur Bestimmung einer rationalen Näherung der Quadratwurzel einer Zahl. Es ist ein Spezialfall des Newton-Verfahrens.

Die Iterationsvorschrift lautet:

$x_{n+1}={\frac {x_{n}+{\frac {a}{x_{n}}}}{2}}$.

Hierbei steht $a$ für die Zahl, deren Quadratwurzel bestimmt werden soll. Der Startwert $x_{0}$ der Iteration kann, solange er nicht gleich Null ist, beliebig festgesetzt werden, wobei zu beachten ist, dass negative Werte gegen die negative Quadratwurzel konvergieren.

Triviales Beispiel

Im Folgenden ein triviales Beispiel für die Wurzel aus 9 und die Annäherung nach vier Berechnungsschritten an den wahren Wert ${\sqrt {9}}=3$

$a=9$ und $x_{0}=1$

$x_{1}={\frac {1+{\frac {9}{1}}}{2}}={\frac {10}{2}}=5$

$x_{2}={\frac {5+{\frac {9}{5}}}{2}}={\frac {\frac {34}{5}}{2}}={\frac {34}{10}}=3{,}4$

$x_{3}={\frac {{\frac {34}{10}}+{\frac {9}{\frac {34}{10}}}}{2}}={\frac {{\frac {34}{10}}+{\frac {90}{34}}}{2}}={\frac {257}{85}}\approx 3{,}0235$

$x_{4}={\frac {{\frac {257}{85}}+{\frac {9}{\frac {257}{85}}}}{2}}={\frac {{\frac {257}{85}}+{\frac {765}{257}}}{2}}={\frac {65537}{21845}}\approx 3{,}00009.$

Geometrische Veranschaulichung des Heron-Verfahrens

Der Flächeninhalt eines Quadrates kann über das Quadrat der Länge seiner Seiten berechnet werden. Die Bestimmung der Quadratwurzel einer gegebenen Zahl $a$ kann also geometrisch gedeutet werden als Bestimmung der Seitenlänge (genauer: als rationale Näherung zu ${\sqrt {a}}$) eines Quadrates mit dem Flächeninhalt $a$.

Die Idee ist nun, von einem Rechteck des Flächeninhaltes $a$ auszugehen und die Seitenlängen einem Quadrat immer weiter anzunähern.

Dazu wird ein Startwert gewählt, im obigen Fall gilt $a=9$ und als Startwert wurde $x_{0}=1$ gewählt. Geometrisch bedeutet dieses, dass von einem Rechteck der Seitenlänge $x_{0}=1$ ausgegangen wird. Die andere Seitenlänge dieses Rechtecks ergibt sich aus dem vorgegebenen Flächeninhalt: $y_{0}={\frac {9}{1}}=9$

Bei der Betrachtung ist unmittelbar ersichtlich, dass es eine geeignetere Näherung an ein Quadrat gibt, denn die eine Seitenlänge $x_{0}=1$ ist zu klein, die andere mit $y_{0}=9$ zu groß.

Um eine verbesserte Annäherung an die Länge einer Quadratseite zu erhalten, kann das arithmetische Mittel der Seitenlängen $x_{0}$ und $y_{0}$ dienen (Hier gibt es eine ganze Reihe von Möglichkeiten, das Verfahren zu verfeinern.):

$x_{1}={\frac {x_{0}+y_{0}}{2}}={\frac {9+1}{2}}=5$

Die Länge der zweiten Seite dieses neuen Näherungs-Rechtecks ergibt sich wieder durch den vorgegebenen Flächeninhalt $a$:

$y_{1}={\frac {a}{x_{1}}}={\frac {9}{5}}=1{,}8$

Die Werte $x_{1}=5$ und $y_{1}=1{,}8$ sind geometrisch gedeutet die Seitenlängen eines zweiten Näherungs-Rechtecks. Dieses und die folgenden Rechtecke lassen sich nun weiter verbessern durch erneute Bildung des arithmetischen Mittelwertes als verbesserte Näherung an die Seitenlänge eines Quadrates mit der Seitenlänge ${\sqrt {a}}$.

Konvergenz

Das Verfahren konvergiert relativ rasch. Man erhält innerhalb weniger Schritte eine sehr gute Näherung. Da es sich aus dem Newtonschen Näherungsverfahren ableiten lässt, ist die Konvergenzordnung 2.
Es gelten:

${\sqrt {a}}\leq x_{n+1}\leq x_{n}\quad (n=1,2...)$

und

$\lim _{n\to \infty }x_{n}={\sqrt {a}}$

Fehler

Für den Fehler der Heron-Folge $(x_{n})_{n\geq 1}$ gilt:

${\frac {a}{x_{n}}}\leq {\sqrt {a}}\leq x_{n}$ (Einschließung),

sowie

$x_{n}-{\sqrt {a}}\leq {\frac {1}{2{\sqrt {a}}}}\left(x_{n-1}-{\sqrt {a}}\right)^{2}$ (quadratische Konvergenz)

Implementierung in Software

Das Verfahren eignet sich besonders gut zur Implementierung in Software, wird heute angesichts der breiten Verfügbarkeit numerischer Prozessorhardware aber nur noch selten benötigt. Man kommt dabei nur mit Grundrechenarten aus, s. o.

Wenn dazu noch eine Gleitkommadarstellung mit einem Zweier-Exponenten benutzt wird, wird der Ansatz relativ einfach, als Beispiel wird die Wurzel aus 5 betrachtet und der relative Fehler zum Endwert (also abs((x_i - x) / x) ) verfolgt:

• Zunächst spaltet man von diesem Zweier-Exponenten eine gerade Anzahl ab, so dass als Exponent entweder eine -1 oder 0 übrigbleibt, die Zahl also auf das Intervall { ½ , 2 } normalisiert wird. In diesem Intervall ist die Wurzelfunktion eine nur schwach gekrümmte Kurve, lässt sich also numerisch gut behandeln. Beispiel: ${\sqrt {5}}={\sqrt {4\cdot 1,25}}=2\cdot {\sqrt {1,25}}\approx 2\cdot 1,118034=2,236068$, es wird also vorerst nur noch a=1,25 mit dem Ziel x=1,118034 behandelt.
• Als Startwert für die eigentliche Iteration approximiert man diese Kurve durch eine noch einfachere, die man direkt (ohne Iteration) berechnen kann. Mit dieser Anfangsberechnung ermittelt man also den Startwert, mit dem man in die folgende Iteration eintritt. Man kann diese Kurve mehr oder weniger aufwendig ansetzen, mit den steigend komplizierteren Ansätzen unten lässt sich ggf. ein Iterationsschritt einsparen:
  • eine einfache Konstante (z. B. 1),
    Beispiel: x₀ = 1; rel. Fehler=1,1*10^-1;
  • eine Gerade mit Steigung 1/2 und einer additiven Konstante von 1/2 (als Vereinfachung des nachfolgenden Falls),
    Beispiel: x₀=1/2+1,25/2=1,125; rel. Fehler=6,2*10^-3;
  • eine Gerade mit Steigung 1/2 und einer additiven, optimierten Konstante von $\left(2\cdot {\sqrt[{4}]{2}}-{\sqrt {2}}\right)^{2}/2\approx 0{,}929683/2$ ,
    Beispiel: x₀=0,929683/2+1,25/2=1,089841; rel. Fehler=2,5*10^-2;
  • eine Gerade mit optimierter Steigung und einer additiven Konstante.
• Ausgehend von dem so ermittelten Startwert führt man eine feste Anzahl von Iterationsschritten durch. Die nötige Anzahl, um die gewünschte Genauigkeit zu erreichen, lässt sich dank der obigen Fehlerabschätzung als Worst Case innerhalb des Startintervalls direkt ausrechnen. Bei 32 Bits Mantisse und dem mittleren Startansatz braucht man z. B. nur drei Schritte. Diese fest gewählte Anzahl erspart wesentlich aufwendigere Abfragen auf Erreichung der Genauigkeit. Der Ersatz der genannten Konstanten durch die Zahl 1,0 ändert daran nichts. Auch der noch kompliziertere Ansatz brächte zumindest bei dieser Genauigkeit keine Einsparung eines weiteren Iterationsschritts. Bei höheren Genauigkeitsanforderungen kann das anders aussehen.
  Beispiel mit drei Schritten nach Ansatz 1 (Konstante 1, mit den anderen Ansätzen konvergiert es noch einen Schritt schneller):
  x₁=(x₀+a/x₀)/2=(1+1,25/1)2=1,125; rel. Fehler=6,2*10^-3
  x₂=(x₁+a/x₁)/2=(1,125+1,25/1,125)/2=1,118056; rel. Fehler=2,0*10^-5
  x₃=(x₂+a/x₂)/2=(1,118056+1,25/1,118056)/2=1,118034; rel. Fehler=0
  Man sieht die Wirkung der quadratischen Konvergenz, dass sich nämlich der Fehler von Schritt zu Schritt jeweils quadriert oder die Anzahl gültiger Stellen bzw. der negative Fehlerexponent grob verdoppelt.
• Zum Schluss muss man den Exponenten restaurieren, indem man die Hälfte des im ersten Schritt abgespalteten Werts wieder hinzufügt.
  Beispiel: x =2 * x₃ = 2,236068 .

Eine Umsetzung in einer C-Funktion:

// Quadratwurzel-Funktion
// wofür?...
// was tut es: Quadratwurzel mit Heron-Alg. berechnen
// gpl
// von hørst nach Wikipedia-Artikel 
//
// Anmerkung: Funktioniert nur, wenn der float type des Compilers 32 Bit groß ist.
//            Herausfinden kann man dies z.B. so: printf("%d\n", sizeof(float) * 8);

float Wurzel(float Input) 
{ 

// float number in ihren 32 Bits:                           // siehe Artikel Gleitkommazahl
// b    bbbbbbbb        bbbbbbbbbbbbbbbbbbbbbbb 
// Vorzeichen
//      Exponent+127(Bias)
//                      Mantisse als ganze Zahl
//                      (Mantisse /(2^Exponent) -1)*2^23     

// Variablen

unsigned long int *Maske;                                   // Maske, um die Bits direkt zu manipulieren
unsigned short int Exponent;                                // Exponenten-Zwischenlager
float Wurzel = 1;                                           // Startwert und Endspeicher
short int i;
Maske = (unsigned long int *)&Input;                        // Maske auflegen 

if (*Maske >>31 == 1)                                       // Vorzeichentest
{	*Maske = 0x7fffffff;
       return Input;}
// Bei negativem Input wird abgebrochen und die für "nan" -> "not a number" reservierte Bitfolge zurückgegeben

if ((*Maske &0xffffff) == 0x7fffff) *Maske +=1;             // Extremfall-Korrektur
// Die Iteration versagt bei ungeradem Exponent und einer Mantisse die binär ausschließlich aus Einsen besteht.(overflow)
// Wurzel aus "fast 1" wäre unkorrigiert = 0.75 

Exponent = *Maske>>23;                                      // Exponent abspeichern 
// *Maske:      bBBBBBBBBbbbbbbbbbbbbbbbbbbbbbbb
// *Maske>>23:                         bBBBBBBBB(bbbbbbbbbbbbbbbbbbbbbbb)
// Exponent =                   00000000BBBBBBBB
// Nur die großen "B"-Bits sind der Exponent, die kleinen "b" sind Mantisse und Vorzeichen (s. o.).
// Eine Bitverschiebung um 23 Bits nach rechts schneidet gerade die Mantisse ab.
// Ob das Vorzeichenbit == 0 ist, wurde schon geprüft.
// Der Exponent wurde als "short unsigned int" deklariert und hat darum nur 16 Bits.

*Maske=((*Maske|0xff000000) &0x3fffffff);                   // einen geraden Anteil vom Exponent abspalten
// *Maske:      bBBBBBBBBbbbbbbbbbbbbbbbbbbbbbbb
// |0xff000000: 11111111000000000000000000000000
// &0x3fffffff: 00111111111111111111111111111111
// *Maske=      00111111Bbbbbbbbbbbbbbbbbbbbbbbb
// 
// Die bitweise ODER-Verknüpfung mit Hex:ff000000 setzt die ersten acht Bit auf Eins, und lässt
// das letzte Bit des Exponenten sowie die Mantisse intakt.
// Die bitweise UND-Verknüpfung mit Hex:3fffffff setzt die ersten beiden Bits auf Null.
// Der Exponent in der Eingangszahl "Input" ist jetzt: 0111111B = 126 oder 127.
// Je nach dem verbliebenen Bit und inklusive Bias (s. o.) entspricht das gerade:
// "2^0" oder "2^-1".
// Input ist jetzt auf {0.5,2} normiert 

Maske=(long unsigned *)&Wurzel;                             // Maske aufs Zwischenergebnis auflegen

// Iteration für die normierte Zahl
for(i=0;i<4;i++)                                            // eine feste Anzahl von Iterationsschritten
       Wurzel = ( Wurzel + Input/Wurzel ) /2;  

*Maske=(*Maske &0x7fffff)+(((*Maske>>23) + Exponent)<<22);  // Exponenten restaurieren
// (
// *Maske:       00111111Bbbbbbbbbbbbbbbbbbbbbbbb
// &0x7fffffff:  00000000011111111111111111111111
// (=):          000000000bbbbbbbbbbbbbbbbbbbbbbb
// )
// +
// (
// *Maske>>23:                          00111111B(bbbbbbbbbbbbbbbbbbbbbbb)
// +Exponent:                    00000000BBBBBBBB
// (=):                          0000000BBBBBBBB0  (*)
// <<22: (000000)0BBBBBBBB00000000000000000000000
// )
// ______________________________________________
//  =            0BBBBBBBBbbbbbbbbbbbbbbbbbbbbbbb
// 
// Der erste Term der Addition ist das Zwischenergebnis mit gelöschtem Exponent.
// Für den zweiten Term wird das Zwischenergebnis wieder 23 Bit nach rechts verschoben, so dass nur der 
// Normexponent bleibt, und dieser wird mit dem ursprünglichen Exponenten addiert. Das weicht etwas 
// vom oben beschriebenen Algorithmus ab, warum, steht weiter unten.
// (*) Dieser Term ist immer grade, weil das letzte Bit vom Normexponent und vom Ursprungsexponent
// identisch ist. Daher ist das letzte Bit des Terms immer Null.
// Der zweite Term muss noch halbiert und 23 Bits weiter links addiert werden.
// Eine Division durch Zwei bei binären Ganzzahlen entspricht einer 1-Bit-Verschiebung nach rechts.
// Weil das letzte Bit immer Null ist (additionsneutral), können die Schritte als Bit-Verschiebung um
// (nur) 22 Bits nach links und Addition zum ersten Term zusammengefasst werden, ohne die Mantisse zu verändern.
// Diese leicht andere Restauration des Exponenten ist dem (genialen) Bias geschuldet.
// Der Bias ist ein fester Wert, der zu den Exponenten addiert ist, so dass sie immer positiv sind.
// Bei "float" ist dieser Bias 127.
// Ihn zu entfernen, auf das Vorzeichen des Ursprungsexponenten zu prüfen, um "richtigrum" zu halbieren und
// wieder einzufügen, kann mit folgender Umformung umgangen werden:
//
//neuer Exponent =  Normexponent + gerader,abgespaltener Teilexponent /2
//               =  Normexponent + (Ursprungsexponent - Normexponent)/2
//               = (Normexponent + Normexponent + Ursprungsexponent - Normexponent)/2
//               = (Normexponent + Ursprungsexponent)/2

return Wurzel;
} 

Ohne Kommentare sind das nur elf Zeilen mit Anweisungen, die meist nur Bits herumschubsen. So eine Funktion läuft auch auf kleinsten Geräten, ohne die große Mathebibliothek zu bemühen. Andererseits ist sie ungenauer:

$\max \left(\left|{\bigl (}{\text{Wurzel}}\left(x\right){\bigr )}^{2}-x\right|/x\right)\approx 0{,}0000001192$

und etwa viermal langsamer als die in <math.h> enthaltene sqrt()-Funktion.

Verallgemeinerung des Verfahrens

Dieses Verfahren kann man verallgemeinern, so dass man die k-te Wurzel von a berechnen kann. Je größer k aber ist, umso mehr Schritte werden benötigt, um die Wurzel genau zu berechnen.

Man muss dazu die k-dimensionale Entsprechung der oben genannten geometrischen Herleitung benutzen.

$x_{n+1}={\frac {x_{n}+{\sqrt[{k-1\,}]{\frac {a}{x_{n}}}}}{2}}$

Man muss dazu aber die k-1-te Wurzel berechnen. Dazu kann man wieder die Verallgemeinerung benutzten, bis man bei der Quadratwurzel ist.

Geschichte

Dieses Verfahren ist auch unter dem Namen Heron-Verfahren bekannt. Es wurde nach Heron von Alexandria benannt und entstammt seiner Formelsammlung.