Noch mehr Nachkommastellen

Wer deutlich mehr Nachkommastellen benötigt, kommt nicht drum herum, einen besseren Algorithmus als den auf den vorherigen Seiten hergeleiteten zu benutzen, denn dieser lässt sich nicht wirklich effizient implementieren. Insbesondere das Wurzelziehen benötigt bei zunehmender Genauigkeit deutlich mehr Rechenleistung. Zum Beispiel eignet sich der 1770 von Johann Heinrich Lambert publizierte Kettenbruch, da er leicht mit einfachen arithmetischen Operationen implementiert werden kann:

$\frac{4}{\pi}=1+\frac{1^2}{3+\frac{2^2}{5+\frac{3^2}{7+\frac{4^2}{9+\frac{5^2}{\ddots}}}}}$

Laut Wikipedia bringt jeder Schritt etwa 0,76555 Dezimalstellen. Ich habe diesen Algorithmus in der funktionalen Programmiersprache Haskell implementiert und bin davon ausgegangen, dass jeder Schritt 0,75 Dezimalstellen bringt. Außerdem sollten prinzipiell erstmal zehn Schritte für die "nullte" Dezimalstelle gemacht werden, damit die Stellen auch mit Sicherheit stimmen:

Haskell-Code
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15:	module Pi (approximate) where import ShowRational -- Berechnet Pi mindestens auf n Stellen genau approximate :: Integer -> Rational approximate n = 4 / (lambert 1 steps) where steps = floor (4/3(fromInteger n) + 10) -- Berechnet den Lambertkettenbruch bis ins nte Glied lambert :: Integer -> Integer -> Rational lambert k n \| (k == n) = 1 \| otherwise = a + b / (lambert (k+1) n) where a = fromInteger (2k-1) b = fromInteger (k^2)

Haskell-Code

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:

module Pi (approximate) where
 
import ShowRational
 
-- Berechnet Pi mindestens auf n Stellen genau
approximate :: Integer -> Rational
approximate n = 4 / (lambert 1 steps)
                where steps = floor (4/3*(fromInteger n) + 10)
 
-- Berechnet den Lambertkettenbruch bis ins nte Glied
lambert :: Integer -> Integer -> Rational
lambert k n | (k == n)  = 1
            | otherwise = a + b / (lambert (k+1) n)
            where a = fromInteger (2*k-1)
                  b = fromInteger (k^2)

Folgendes Modul ist nur für die Gleitkomma-Ausgabe rationaler Zahlen (die in Haskell als gekürzter Bruch dargestellt werden) nötig (gerundet wird nicht):

Haskell-Code
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14:	module ShowRational (showRational) where -- Wandelt r mit n Nachkommastellen in eine Zeichenkette um showRational :: Integer -> Rational -> String showRational n r = (show intDigits) ++ "." ++ (showFracs n fracDigits) where intDigits = truncate r fracDigits = 10 * (r-(fromInteger intDigits)) -- Wandelt die ersten n Nachkommastellen in eine Zeichenkette showFracs :: Integer -> Rational -> String showFracs 0 r = "" showFracs n r = (show first) ++ (showFracs (n-1) rest) where first = truncate r rest = 10 * (r - (fromInteger first))

Haskell-Code

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:

module ShowRational (showRational) where
 
-- Wandelt r mit n Nachkommastellen in eine Zeichenkette um
showRational :: Integer -> Rational -> String
showRational n r = (show intDigits) ++ "." ++ (showFracs n fracDigits)
                   where intDigits  = truncate r
                         fracDigits = 10 * (r-(fromInteger intDigits))
 
-- Wandelt die ersten n Nachkommastellen in eine Zeichenkette
showFracs :: Integer -> Rational -> String
showFracs 0 r = ""
showFracs n r = (show first) ++ (showFracs (n-1) rest)
              where first = truncate r
                    rest  = 10 * (r - (fromInteger first))

Zum Beispiel mit dem Aufruf

showRational 100 (approximate 100)

berechnet dieses Programm Pi auf 100 Nachkommastellen genau und gibt das Ergebnis dann auch auf 100 Nachkommastellen genau aus. Auf meinem Laptop (CPU: AMD Turion™ 64 X2 TL-50 @ 1,6 GHz (Dual-Core) mit dem Haskell-Interpreter Hugs, wobei nur einer der beiden Prozessoren genutzt wurde) benötigte dieser Algorithmus ohne die Gleitkomma-Umwandlung (also nur der Aufruf "approximate 1000") knapp 4:30 Minuten, um Pi auf 1000 Nachkommastellen genau zu berechnen: Ergebnis (mit Gleitkomma-Umwandlung).

Zurück: Umsetzung durch ein C++ Programm