Skip to Main Content Skip to Search
Accelerating the pace of engineering and science

 

Cleve’s Corner – Viele Wege führen zu π

Von Cleve Moler

Die Berechnung der Zahl π auf mehrere hundert oder sogar Billionen von Stellen war lange ein probates Mittel, Stresstests an Hardware vorzunehmen, Software zu validieren oder ganz einfach auf sich aufmerksam zu machen. Die MATLAB®-Implementierungen der am häufigsten zur Berechnung von π eingesetzten Algorithmen sind exemplarisch für zwei ganz unterschiedliche Typen von Arithmetik, die in der Symbolic Math Toolbox™ zur Verfügung stehen: exakte rationale Arithmetik und Fließkommaarithmetik mit variabler Genauigkeit.

Im vergangenen August verkündeten die Hobbyinformatiker Alexander Yee und Shigeru Kondo einen neuen Weltrekord mit der Bestimmung von 5 Billionen Stellen von π (Abbildung 1). Ihre Berechnung dauerte 90 Tage auf einem Computer-Eigenbau. Yee absolviert derzeit ein Graduiertenstudium an der University of Illinois. Seine Rechensoftware „y-cruncher“ entstand aus einem Schulprojekt an der Palo Alto High School. Kondo ist Systemingenieur und baut sich zu Hause in Japan seine PCs selbst auf. Der Rechner, den er für dieses Projekt zusammengestellt hat (Abbildung 2), verfügt über zwei Intel® Xeon®-Prozessoren mit insgesamt 12 Kernen, 96 Gigabyte RAM und 20 externe Festplatten mit einer Gesamtkapazität von 32 Terabyte.

Abb. 1: Logarithmische Darstellung der im Wikipedia-Artikel „Chronology of computation of pi“ aufgelisteten Daten. Sie zeigt, wie sich der Weltrekord der Anzahl berechneter Stellen in den letzten 60 Jahren entwickelt hat. Die lineare Anpassung an den Logarithmus verrät, dass hier eine Variante des Mooreschen Gesetzes erfüllt ist. Die Stellenzahl verdoppelt sich im Schnitt alle 22,5 Monate.

Abb. 1: Logarithmische Darstellung der im Wikipedia-Artikel „Chronology of computation of pi“ aufgelisteten Daten. Sie zeigt, wie sich der Weltrekord der Anzahl berechneter Stellen in den letzten 60 Jahren entwickelt hat. Die lineare Anpassung an den Logarithmus verrät, dass hier eine Variante des Mooreschen Gesetzes erfüllt ist. Die Stellenzahl verdoppelt sich im Schnitt alle 22,5 Monate.

Abb. 2: Shigeru Kondos Eigenbau-PC, der derzeit den Weltrekord der berechneten Stellen von π hält.

Abb. 2: Shigeru Kondos Eigenbau-PC, der derzeit den Weltrekord der berechneten Stellen von π hält.

Die Berechnung von 5 Billionen Stellen ist eine gigantische Datenverarbeitungsaufgabe, bei der die benötigte Zeit für den Datentransfer genauso entscheidend ist wie die eigentliche Rechenzeit.

y-cruncher nutzt mehrere verschiedene Formeln zur Berechnung und Verifikation. Sein Herzstück ist eine eindrucksvolle Formel, die 1987 von David und Gregory Chudnovsky entdeckt wurde:

1 π = 12 k = 0 - 1 k 6 k ! 13591409 + 545140134 k 3 k ! k ! 3 640320 3 k + 3 / 2

Die Brüder Chudnovsky, denen zwei interessante Artikel im Magazin The New Yorker sowie ein NOVA-Dokumentarfilm gewidmet wurden, haben ebenfalls ihre eigenen Computer in ihrer Wohnung in Manhattan zusammengebaut. Der Massenparallelrechner, den sie in den 90er Jahren aus handelsüblichen Teilen zusammengestellt haben, galt damals als Supercomputer.

Hochexakte Arithmetik in der Symbolic Math Toolbox

Die Symbolic Math Toolbox bietet zwei Arten von hochexakter Arithmetik: sym und vpa.

Die Funktion sym initiiert die exakte rationale Arithmetik: Arithmetische Größen werden durch Quotienten sowie durch Wurzeln großer Ganzzahlen dargestellt. Die Integer-Länge wird je nach Bedarf vergrößert und ist ausschließlich durch die Rechenzeit und den benötigten Speicher begrenzt. Quotienten und Wurzeln werden nur berechnet, wenn das Ergebnis eine exakte Ganzzahl ist.

Die Funktion vpa dagegen ruft die Fließkommaarithmetik mit variabler Genauigkeit auf: Hier werden arithmetische Größen durch Dezimalbrüche mit einer bestimmten Stellenzahl sowie eine angehängte Zehnerpotenz dargestellt.

Arithmetische Operationen, auch Divisionen und Wurzeln, können Rundungsfehler aufweisen, deren Größe durch die gewählte Genauigkeit bestimmt ist. So erzeugt etwa die MATLAB-Anweisung

p = rat(pi)

drei Terme der Näherung von π durch einen Kettenbruch

p = 3 + 1/(7 + 1/16)

Die Anweisung

sym(p)

liefert danach den rationalen Ausdruck

355/113

der eine reizvolle Alternative zum sonst üblichen Bruch 22/7 darstellt. Dagegen ergibt die Anweisung

vpa(p,8)

den achtstelligen Fließkommawert

3.1415929

Der Chudnovsky-Algorithmus

Die vorstehenden MATLAB-Programme implementieren nur drei aus hunderten möglicher Algorithmen zur Berechnung von π. Der Chudnovsky-Algorithmus (Abbildung 3) bedient sich bis auf den letzten Schritt der exakten rationalen Arithmetik.

function P = chud_pi(d)
% CHUD_PI Chudnovsky algorithm for pi.
% chud_pi(d) produces d decimal digits.

k = sym(0);
s = sym(0);
sig = sym(1);
n = ceil(d/14);
for j = 1:n
    s = s + sig * prod(3*k+1:6*k)/prod(1:k)^3 * ...
    (13591409+545140134*k) / 640320^(3*k+3/2);
    k = k+1;
    sig = -sig;
end
S = 1/(12*s);
P = vpa(S,d);

Abb. 3: MATLAB-Implementierung des Chudnovsky-Algorithmus. sym initiiert darin die exakte rationale Arithmetik.

Jeder Ausdruck der Chudnovsky-Reihe liefert etwa 14 Dezimalstellen, d. h. um d Stellen zu erhalten, benötigt man [d/14] Ausdrücke. Die Anweisung

chud_pi(42)

erzeugt beispielsweise durch exakte rationale Berechnung

3405449678154416971853498155008000000 640320 867407410133324147761288805130794983129

Der letzte Schritt erzeugt dann mithilfe von Fließkommaarithmetik mit variabler Genauigkeit das 42-stellige Ergebnis

3.14159265358979323846264338327950288419717

Eine interessante Besonderheit in den Ziffern von π offenbart sich durch

chud_pi(775)

oder

vpa(pi,775)

Beide erzeugen 775 Stellen, die jeweils anfangen und enden mit:

3.141592653589793238 ...... 707211349999998372978

Wir erkennen hier nicht nur, dass chud_pi korrekt arbeitet, sondern auch, dass um die Position 765 in der Dezimalentwicklung von π sechsmal hintereinander die Ziffer 9 erscheint (Abbildung 4).

Abb. 4: 10.000 Stellen von π in einer MATLAB-Grafik. Erkennen Sie die sechs aufeinanderfolgenden 9er in der achten Reihe?

Abb. 4: 10.000 Stellen von π in einer MATLAB-Grafik. Erkennen Sie die sechs aufeinanderfolgenden 9er in der achten Reihe?

Die Berechnung von 5000 Stellen mit chud_pi erfordert 358 Ausdrücke und dauert auf meinem Laptop ungefähr 20 Sekunden. Der dabei erzeugte symbolische Ausdruck enthält 12.425 Zeichen.

Der Algorithmus des algebraisch-geometrischen Mittels

Die Chudnovsky-Formel ist eine Potenzreihe: Jeder neue Ausdruck der Partialsumme vergrößert die Stellenzahl um einen bestimmten Betrag. Der Algorithmus für das algebraisch-geometrische Mittel funktioniert völlig anders; er ist quadratisch konvergent: Jede neue Iteration verdoppelt die Zahl der korrekten Stellen. Dieser Algorithmus hat eine lange Geschichte, die bis auf Gauss und Legendre zurückgeht. Dass er sich auch zur Berechnung von π eignet, haben 1975 Richard Brent und Eugene Salamin unabhängig voneinander entdeckt.

Das arithmetische Mittel zweier Zahlen, (a+b)/2, ist immer größer als ihr geometrisches Mittel, a b . Ihr arithmetisch-geometrisches Mittel, agm(a,b), wird durch wiederholte Ermittlung arithmetischer und geometrischer Mittelwerte berechnet. Ausgehend von a0 = a und b0 = b, iteriert man

a n + 1 = a n + b n 2 b n + 1 = a n b n

bis an und bn mit gewünschter Genauigkeit übereinstimmen. Das arithmetische Mittel von 1 und 9 ist beispielsweise 5, ihr geometrisches Mittel 3 und das agm ist 3,9362.

Die mächtigen Eigenschaften des agm liegen nicht nur im Endergebnis, sondern auch in den errechneten Zwischenergebnissen. Die Berechnung von π ist nur ein Beispiel dafür. Man berechnet dazu

M = agm ( 1 , 1 2 )

Während der Iteration überwacht man die Änderungen

c n = a n + 1 - a n

Schließlich sei

S=14 n = 0 2 n c n2

Hieraus folgt

π = M 2 S

Unsere MATLAB-Funktion für den agm-Algorithmus (Abbildung 5) arbeitet von Anfang an mit Fließkommaarithmetik mit variabler Genauigkeit. Würden wir symbolische rationale Arithmetik verwenden, dann erhielten wir eine verschachtelte Reihe von Quadratwurzeln. Nach nur zwei Iterationen erhielten wir

2 2 2 + 2 8 + 1 4 2 1 4 - 2 4 - 1 2 2 - 2 2 8 - 2 2 2 + 1 4 2

Der Dezimalwert dieses Ungetüms ist 3,14168 und noch keine besonders gute Annäherung an π. Die Komplexität verdoppelt sich mit jeder weiteren Iteration.

function P = agm _ pi(d)
% AGM _ PI Arithmetic-geometric mean for pi.
% agm _ pi(d) produces d decimal digits.
digits(d)
a = vpa(1,d);
b = 1/sqrt(vpa(2,d));
s = 1/vpa(4,d);
p = 1;
n = ceil(log2(d));
for k = 1:n
    c = (a+b)/2;
    b = sqrt(a*b);
    s = s - p*(c-a)^2;
    p = 2*p;
    a = c;
end
P = a^2/s;

Abb. 5: MATLAB-Implementierung des Algorithmus für das arithmetisch-geometrische Mittel. digits und vpa initiieren die Fließkommaarithmetik mit variabler Genauigkeit.

Solche exakten symbolischen Ausdrücke sind rechnerisch unhandlich und ineffizient. Bei vpa dagegen erkennt man sofort die quadratische Konvergenz an den aufeinanderfolgenden Iterationen:

    3.2
    3.142
    3.1415927
    3.141592653589793
    3.1415926535897932384626433832795
    3.1415926535897932…79502884197169399375105820974944592

Der sechste Eintrag in der Liste hat bereits 64 korrekte Stellen, der zehnte 1024.

Durch seine quadratische Konvergenz ist agm sehr schnell. Um auf 100.000 Stellen zu kommen, benötigt mein Laptop nur 17 Schritte und ca. 2,5 Sekunden.

Die BBP-Formel

Yee hat seine Berechnungen stichprobenartig mit der BBP-Formel überprüft. Diese nach David Bailey, Jonathan Borwein und Simon Plouffe benannte und von Plouffe im Jahr 1995 entdeckte Formel ist eine Potenzreihe, die auf den inversen Potenzen von 16 fußt:

π = k = 0 4 8 k + 1 - 2 8 k + 4 - 1 8 k + 5 - 1 8 k + 6 1 16 k

Bemerkenswert an der BBP-Formel ist, dass man mit ihr mehrere aufeinanderfolgende Hexadezimalstellen in der Hexadezimalerweiterung von π berechnen kann, ohne die vorherigen Stellen zu kennen und ohne die Genauigkeit mithilfe zusätzlicher arithmetischer Methoden festzulegen. Eine Formel für die Hexadezimalstellen ab der Position d erhält man ganz einfach, indem man die BBP-Formel mit 16d multipliziert und nur die Nachkommastellen des Ergebnisses betrachtet.

Ich kann nicht widerstehen, dieses Corner mit der Bemerkung zu schließen, dass uns BBP immer nur ein Stückchen des Gesamtbildes von π zeigt, sozusagen „just a piece of π“ (pie: Kuchenstück).

Unser BBP-Programm ist zu lang, um hier aufgeführt zu werden, aber die drei Programme chudovsky_pi, agm_pi, und bbp_pi können auf MATLAB Central heruntergeladen werden.

Veröffentlicht 2011
91964v00

Literatur

Eingesetzte Produkte

Weitere Informationen