Heron-Verfahren

Das Heron-Verfahren (seltener auch Babylonisches Wurzelziehen genannt) 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, 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 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, d. h. es erzeugt mit wenigen Schritten eine 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, da nur Grundrechenarten benötigt werden, s. o. Es wird heute angesichts der breiten Verfügbarkeit numerischer Prozessorhardware aber nur noch selten benötigt.

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 wird von diesem Zweier-Exponenten eine gerade Anzahl abgespaltet, 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 sich direkt (ohne Iteration) berechnen lässt. Mit dieser Anfangsberechnung wird der Startwert ermittelt, mit dem die folgende Iteration begonnen wird. 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 (beipielsweise 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 beispielsweise 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 der Fehler von Schritt zu Schritt jeweils quadriert oder die Anzahl gültiger Stellen bzw. der negative Fehlerexponent grob verdoppelt.
Zum Schluss wird der Exponent restauriert, 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 nach dem Heron-Alg.
 *
 * gpl
 * von hørst nach Wikipedia-Artikel 
 *
 * Voraussetzung für dieses Beispiel ist ein Rechner auf Little-Endian-Basis
 * und mit IEEE754-Gleitkommazahlen.
 *
 * Muss mit -fno-strict-aliasing kompiliert werden.
 */

#include <stdio.h>
#include <stdint.h>

/*
 * IEEE754-32bit-Gleitkommazahl (Voraussetzung):
 *
 * b    bbbbbbbb        bbbbbbbbbbbbbbbbbbbbbbb 
 * Vorzeichen
 *      Exponent+127(Bias)
 *                      Mantisse als ganze Zahl
 *                      (Mantisse / (2^Exponent) - 1) * 2^23
 */

static float wurzel(float input)
{
	unsigned short exponent; /* Exponenten-Zwischenlager */
	uint32_t *maske;         /* Maske, um die Bits direkt zu manipulieren */
	float wurzel = 1;        /* Startwert und Endspeicher */
	short i;

	/* Maske auflegen */
	maske = (uint32_t *)&input;

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

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

	/*
	 * 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" deklariert
	 * und hat darum mindestens 16 Bits.
	 */
	exponent = *maske >> 23;


	/*
	 * 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 0x3FFFFFFF 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 = ((*maske | 0xFF000000) & 0x3FFFFFFF);

	/* Maske aufs Zwischenergebnis auflegen */
	maske = (uint32_t *)&wurzel;

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


	/* 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 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
	 */
	*maske = (*maske &0x7fffff) + (((*maske >> 23) + exponent) << 22);

	return wurzel;
}

int main(void)
{
	/*
	 * Because the original author decided to write an
	 * endian-dependent example. (Operating on a float using
	 * (unsigned long *)/(uint32_t *))
	 */
	float neg = -0.0;
	if (*(unsigned char *)&neg != 0)
		printf("Example will not correctly work on this processor\n");

	printf("%f\n", wurzel(2));
	return 0;
}

Ohne Kommentare sind das nur elf Zeilen mit Anweisungen, die meist nur Bits herumschubsen. Eine solche 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 lässt sich verallgemeinern, so dass die $k$ -te Wurzel von $a$ berechnet wird. Je größer $k$ ist, desto mehr Schritte werden benötigt, um die Wurzel genau zu berechnen.

Dazu wird die k-dimensionale Entsprechung der oben genannten geometrischen Herleitung benutzt.

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

Dazu wird die $k-1$ -te Wurzel berechnet, wozu wieder die Verallgemeinerung benutzt wird, bis die Quadratwurzel erreicht ist.

Geschichte

Das Verfahren wurde nach Heron von Alexandria benannt und entstammt seiner Formelsammlung.