Chudnovsky-Algorithmus

Der Chudnovsky-Algorithmus ist eine von den Chudnovsky-Brüdern im Jahre 1988 entwickelte iterative Methode zur Berechnung beliebig vieler Nachkommastellen der Kreiszahl π. Jede Iteration liefert durchschnittlich 14,81 weitere Dezimalstellen. Der Algorithmus basiert auf der Konvergenz einer verallgemeinerten hypergeometrischen Reihe:

${\frac {1}{\pi }}=12\sum _{k=0}^{\infty }{\frac {(-1)^{k}(6k)!(545140134k+13591409)}{(3k)!(k!)^{3}640320^{3k+3/2}}}.\!$

Dieser Algorithmus wurde seitdem für alle Weltrekordberechnungen eingesetzt: 2,7 Billionen Stellen von π im Dezember 2009; 5 Billionen Stellen von π im August 2010; 10 Billionen Stellen von π im Oktober 2011; 12 Billionen Stellen von π im Dezember 2013; 31,4 Billionen Stellen von π im Januar 2019; 50 Billionen Stellen von π im Januar 2020, über 62,8 Billionen Stellen im August 2021 und 100 Billionen im Juni 2022.

Entwicklung

Heegner-Punkte können dabei helfen, sehr schnell konvergente Reihen zu finden, die gegen die Kreiszahl $\pi$ konvergieren. Vorläufer solcher Reihentypen waren schon von Srinivasa Ramanujan zu Beginn des 20. Jahrhunderts entdeckt worden. Die Brüder David und Gregory Chudnovsky nutzten schließlich die Punkte $\tau ={\tfrac {1+i{\sqrt {n}}}{2}}$ mit natürlichen Zahlen $n$, um die Arbeiten von Ramanujan weiterzuführen. Dabei fanden sie eine für $j(\tau ):=1728J(\tau )$ und all diese Heegner-Punkte gültige Reihenidentität

$\sum _{k=0}^{\infty }\left({\frac {1}{6}}(1-s_{2}(\tau ))+k\right){\frac {(6k)!}{(3k)!k!^{3}}}{\frac {1}{j(\tau )^{k}}}={\frac {\sqrt {-J(\tau )}}{\pi }}{\frac {1}{\sqrt {n(1-J(\tau ))}}},$

die den durch Eisensteinreihen definierten Term

$s_{2}(\tau )={\frac {E_{4}(\tau )}{E_{6}(\tau )}}\left(E_{2}(\tau )-{\frac {3}{\pi \mathrm {Im} (\tau )}}\right)$

beinhaltet. Dabei bezeichnet $k!$ die Fakultät von $k$. Daraus konnte nach Einsetzen des Heegner-Punkts $\tau ={\tfrac {1+i{\sqrt {163}}}{2}}$ der Chudnovsky-Algorithmus entwickelt werden, mit Hilfe dessen die Kreiszahl $\pi$ extrem schnell auf viele Nachkommastellen berechnet werden kann. Er nutzt aus, dass der Wert $j\left({\tfrac {1+i{\sqrt {163}}}{2}}\right)$ ganzzahlig ist. Über die Methoden, wie man $j(\tau )$ allgemein berechnet, kann man bereits diese und weitere Kuriositäten beobachten. Man weiß wegen der Fourier-Entwicklung $j(\tau )=e^{-2\pi i\tau }+744+196\,884e^{2\pi i\tau }+\dotsb$, dass für Werte $\tau$ mit größerem Imaginärteil die Zahl $j(\tau )$ bereits sehr nahe an $e^{-2\pi i\tau }+744$ liegt. In der Tat findet man

$e^{\pi {\sqrt {163}}}=262\,537\,412\,640\,768\,743{,}\,999\,999\,999\,999\,250\ldots .$

Ein ausführlicher Beweis dieser Formel findet sich hier:

Diese ist ähnlich der Ramanujan-Formel zur Ermittlung von π, und ist ein Beispiel der Ramanujan-Sato-Reihen.

Implementierung des Algorithmus

$\pi$ kann mit beliebiger Genauigkeit durch Implementierung des oben genannten Algorithmus in eine dafür geeignete Programmierungsumgebung berechnet werden. Im Folgenden ein Beispiel mit MATLAB:

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