Autor: Andre Adrian
Version: 12.Jun.2011
Die Quadratwurzel von 4 ist 2. Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628... Wie berechnet ein Taschenrechner oder ein Computer die Quadratwurzel?
Heron von Alexandria hat im ersten Jahrhundert nach Christi das Heron-Verfahren, auch babylonisches Wurzelziehen genannt, dokumentiert. Mit der Intervallschachtelung (Bisektion) und dem schriftlichen Wurzelziehen (stellenweises Wurzelziehen) gibt es zwei weitere Näherungsverfahren. Für Gleitkommazahlen und ganze Zahlen werden die Näherungsverfahren als C Quelltext angegeben. Die Berechnung von Startwerten wird behandelt. Gute Startwerte sind wichtig.
double sqrt_heron(double a)
{
double x = a; // Startwert = a
for (int i = 1; i <= 5; ++i) {
x = (x + a / x) / 2;
}
return x;
}
Die C-Funktion sqrt_heron() berechnet die Quadratwurzel x für den Übergabeparameter a. Als Anfangswert für das Heron-Verfahren wird a verwendet.
sqrt_heron(545454545)
Start x = 545454545.000000
i = 1 x = 272727273.000000
i = 2 x = 136363637.500000
i = 3 x = 68181820.750000
i = 4 x = 34090914.375000
i = 5 x = 17045465.187499
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
Der Startwert x = a ist schlecht.
double sqrt_heron1(double a) {
int exponent, i;
double fraction = frexp(a, &exponent);
exponent = exponent / 2; /* Startwert ist a mit Exponent halbiert */
double x = ldexp(fraction, exponent);
for (i = 1; i <= 5; ++i) {
x = (x + a / x) / 2;
}
return x;
}
Mit einem guten Startwert kann die Anzahl der nötigen Iterationen reduziert werden. Natürlich sollte ein guter Startwert einfach zu berechnen sein. Bei Gleitkommazahlen bietet es sich an den Exponenten zu halbieren um einen Startwert zu erhalten. Die Zahl 16.0 wird als 1.0 * 2 ^ 4 kodiert. Dabei ist die Mantisse 1.0 und der Exponent ist 4. Wird der Exponent 4 halbiert ergibt sich 1.0 * 2 ^ 2 oder 4.0.
Die C Bibliothek Funktion frexp() zerlegt ein Fließkomma-Zahl in die Komponenten Mantisse (fraction) und Exponent. Die Funktion ldexp() fügt Mantisse und Exponent wieder zu einer Fließkommazahl zusammen. Die Funktionen frexp() und ldexp() sind definiert in math.h.
sqrt_heron1(545454545)
Start x = 16645.951691
i = 1 x = 24706.975845
i = 2 x = 23391.960385
i = 3 x = 23354.997565
i = 4 x = 23354.968315
i = 5 x = 23354.968315
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
Der Startwert x = a mit Exponent halbiert ist gut.
double sqrt_bisection(double a)
{
double u = 0; // unterer Startwert
double o = a; // oberer Startwert
double x;
for (int i = 1; i <= 30; ++i) {
x = 0.5 * (u + o);
if (x*x > a) {
o = x;
} else {
u = x;
}
}
return x;
}
Die Intervallschachtelung oder Bisektion benötigt keine Division. Die Konvergenzgeschwindigkeit ist deutlich schlechter als beim Heron-Verfahren.
sqrt_bisection(545454545)
Start [ 0.000000 545454545.000000]
i = 1 x = 272727272.500000 [ 0.000000 272727272.500000]
i = 2 x = 136363636.250000 [ 0.000000 136363636.250000]
i = 3 x = 68181818.125000 [ 0.000000 68181818.125000]
i = 4 x = 34090909.062500 [ 0.000000 34090909.062500]
i = 5 x = 17045454.531250 [ 0.000000 17045454.531250]
i = 6 x = 8522727.265625 [ 0.000000 8522727.265625]
i = 7 x = 4261363.632813 [ 0.000000 4261363.632813]
i = 8 x = 2130681.816406 [ 0.000000 2130681.816406]
i = 9 x = 1065340.908203 [ 0.000000 1065340.908203]
i = 10 x = 532670.454102 [ 0.000000 532670.454102]
i = 11 x = 266335.227051 [ 0.000000 266335.227051]
i = 12 x = 133167.613525 [ 0.000000 133167.613525]
i = 13 x = 66583.806763 [ 0.000000 66583.806763]
i = 14 x = 33291.903381 [ 0.000000 33291.903381]
i = 15 x = 16645.951691 [ 16645.951691 33291.903381]
i = 16 x = 24968.927536 [ 16645.951691 24968.927536]
i = 17 x = 20807.439613 [ 20807.439613 24968.927536]
i = 18 x = 22888.183575 [ 22888.183575 24968.927536]
i = 19 x = 23928.555555 [ 22888.183575 23928.555555]
i = 20 x = 23408.369565 [ 22888.183575 23408.369565]
i = 21 x = 23148.276570 [ 23148.276570 23408.369565]
i = 22 x = 23278.323067 [ 23278.323067 23408.369565]
i = 23 x = 23343.346316 [ 23343.346316 23408.369565]
i = 24 x = 23375.857941 [ 23343.346316 23375.857941]
i = 25 x = 23359.602128 [ 23343.346316 23359.602128]
i = 26 x = 23351.474222 [ 23351.474222 23359.602128]
i = 27 x = 23355.538175 [ 23351.474222 23355.538175]
i = 28 x = 23353.506199 [ 23353.506199 23355.538175]
i = 29 x = 23354.522187 [ 23354.522187 23355.538175]
i = 30 x = 23355.030181 [ 23354.522187 23355.030181]
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
Die Startwerte u = 0, o = a sind brauchbar.
double sqrt_bisection1(double a) {
int exponent, i;
double fraction = frexp(a, &exponent);
/* unterer Startwert ist Exponent halbiert - 1 */
double u = ldexp(fraction, exponent / 2 - 1);
/* oberer Startwert ist Exponent halbiert + 1 */
double o = ldexp(fraction, exponent / 2 + 1);
double x;
for (i = 1; i <= 30; ++i) {
x = 0.5 * (u + o);
if (x * x > a) {
o = x;
} else {
u = x;
}
printf("%d %f\n", i, x);
}
return x;
}
Mit der Exponent halbieren Methode können der untere und obere Startwert für die Intervallschachtelung bestimmt werden.
sqrt_bisection1(545454545)
Start [ 8322.975845 33291.903381]
i = 1 x = 20807.439613 [ 20807.439613 33291.903381]
i = 2 x = 27049.671497 [ 20807.439613 27049.671497]
i = 3 x = 23928.555555 [ 20807.439613 23928.555555]
i = 4 x = 22367.997584 [ 22367.997584 23928.555555]
i = 5 x = 23148.276570 [ 23148.276570 23928.555555]
i = 6 x = 23538.416063 [ 23148.276570 23538.416063]
i = 7 x = 23343.346316 [ 23343.346316 23538.416063]
i = 8 x = 23440.881189 [ 23343.346316 23440.881189]
i = 9 x = 23392.113753 [ 23343.346316 23392.113753]
i = 10 x = 23367.730035 [ 23343.346316 23367.730035]
i = 11 x = 23355.538175 [ 23343.346316 23355.538175]
i = 12 x = 23349.442246 [ 23349.442246 23355.538175]
i = 13 x = 23352.490211 [ 23352.490211 23355.538175]
i = 14 x = 23354.014193 [ 23354.014193 23355.538175]
i = 15 x = 23354.776184 [ 23354.776184 23355.538175]
i = 16 x = 23355.157180 [ 23354.776184 23355.157180]
i = 17 x = 23354.966682 [ 23354.966682 23355.157180]
i = 18 x = 23355.061931 [ 23354.966682 23355.061931]
i = 19 x = 23355.014306 [ 23354.966682 23355.014306]
i = 20 x = 23354.990494 [ 23354.966682 23354.990494]
i = 21 x = 23354.978588 [ 23354.966682 23354.978588]
i = 22 x = 23354.972635 [ 23354.966682 23354.972635]
i = 23 x = 23354.969658 [ 23354.966682 23354.969658]
i = 24 x = 23354.968170 [ 23354.968170 23354.969658]
i = 25 x = 23354.968914 [ 23354.968170 23354.968914]
i = 26 x = 23354.968542 [ 23354.968170 23354.968542]
i = 27 x = 23354.968356 [ 23354.968170 23354.968356]
i = 28 x = 23354.968263 [ 23354.968263 23354.968356]
i = 29 x = 23354.968310 [ 23354.968310 23354.968356]
i = 30 x = 23354.968333 [ 23354.968310 23354.968333]
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
Die Startwerte u = a mit Exponent halbiert -1 , o = a mit Exponent halbiert +1 sind gut.
Die Rechenwerke im Computer arbeiten intern mit ganzen Zahlen. Die Mantisse einer Fließkomma-Zahl ist eine Festkomma-Zahl, im gewissen Sinne ein Integer.
unsigned int isqrt_heron(unsigned int a)
{
unsigned int b = 1 << 30; // Vergleichswert ist 2 ^ 30
unsigned int x = 1 << 15; // Startwert ist Quadratwurzel von b
// Startwert bestimmen, Vergleichswert b kleiner als a machen
while (b > a) {
b >>= 2; // Vergleichswert um 4 reduzieren
x >>= 1; // Startwert um sqrt(4) reduzieren
}
for (int i = 1; i <= 4; ++i) {
x = (x + a / x) / 2;
}
return x;
}
Bei Integer (ganzen Zahlen) wird ein guter Startwert für das Heron-Verfahren mit Hilfe eines Vergleichswertes bestimmt. In einer while() Schleife wird der Vergleichswert um den Faktor 4 reduziert, der Startwert um den Faktor 2.
isqrt_heron(545454545)
Start x = 16384
i = 1 x = 24837
i = 2 x = 23399
i = 3 x = 23355
i = 4 x = 23354
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
unsigned int isqrt_bisection(unsigned int a)
{
unsigned int b = 1 << 30; // Vergleichswert ist 2 ^30
unsigned int u = 1 << 15; // unterer Startwert ist Wurzel von Vergleichswert
// unteren Startwert bestimmen, Vergleichswert b kleiner als a machen
while (b > a) {
b >>= 2;
u >>= 1;
}
unsigned int o = 2 * u; // oberer Startwert
unsigned int x;
for (int i = 1; i <= 15; ++i) {
x = (u + o) >> 1;
if (x*x > a) {
o = x;
} else {
u = x;
}
}
return x;
}
isqrt_bisection(545454545)
Start [16384 32768]
i = 1 x = 24576 [16384 24576]
i = 2 x = 20480 [20480 24576]
i = 3 x = 22528 [22528 24576]
i = 4 x = 23552 [22528 23552]
i = 5 x = 23040 [23040 23552]
i = 6 x = 23296 [23296 23552]
i = 7 x = 23424 [23296 23424]
i = 8 x = 23360 [23296 23360]
i = 9 x = 23328 [23328 23360]
i = 10 x = 23344 [23344 23360]
i = 11 x = 23352 [23352 23360]
i = 12 x = 23356 [23352 23356]
i = 13 x = 23354 [23354 23356]
i = 14 x = 23355 [23354 23355]
i = 15 x = 23354 [23354 23355]
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
Der folgende Quelltext basiert auf den Quelltext im Wikipedia Artikel Methods of computing square roots
unsigned int isqrt(unsigned int a) {
unsigned int b = 1 << 30; // Vergleichswert beginnt mit 2 ^30
unsigned int x = 0;
for (int i = 1; i <= 16; ++i) {
if (a >= x + b) {
a -= (x + b);
x += (b + b);
}
x >>= 1;
b >>= 2;
}
return x;
}
Das schriftliche Wurzelziehen ist im Binärsystem ohne Multiplikation möglich. Die Konvergenzgeschwindigkeit ist wie bei der Intervallschachtelung.
isqrt(545454545)
x = 0 a = 545454545 b = 268435456
x = 268435456 a = 277019089 b = 67108864
x = 134217728 a = 277019089 b = 16777216
x = 83886080 a = 126024145 b = 4194304
x = 46137344 a = 37943761 b = 1048576
x = 23068672 a = 37943761 b = 262144
x = 11796480 a = 14612945 b = 65536
x = 5963776 a = 2750929 b = 16384
x = 2981888 a = 2750929 b = 4096
x = 1490944 a = 2750929 b = 1024
x = 746496 a = 1258961 b = 256
x = 373504 a = 512209 b = 64
x = 186816 a = 138641 b = 16
x = 93408 a = 138641 b = 4
x = 46708 a = 45229 b = 1
x = 23354 a = 45229 b = 0
Die Quadratwurzel von 545454545 ist 23354,968315114452323438891821628...
Neben Intel mit dem 8086 und Motorola mit dem 68000 hat sich National Semiconductor mit dem NS32016 um den 16-bit/32-bit Mikroprozessor Markt bemüht. Wie bekannt hat Intel mit dem 8088 zuerst die Gunst von IBM gewonnen und später einen Großteil des Mikroprozessor-Marktes.
Der Floating Point Coprozessor NS32081 beherrscht nur die vier Grundrechenarten in Gleitkomma-Arithmetrik. Vor vielen Jahren durfte ich das Gespann NS32016/NS32081 programmieren.
Eine Frage seinerzeit war: was ist die beste Quadratwurzel Routine für diese FPU?