Source: https://www.unet.univie.ac.at/~a8727063/Science/BBP/
Timestamp: 2018-11-20 01:02:04+00:00

Document:
Konkretes Beispiel für den BBP-Algorithmus (Bailey Borwein Plouffe) für Pi mit n = 5, n=1.000 (k=20) und n=(10^9-1)
[http://kurzurl.net/bbp]
Konkretes Beispiel für den BBP-Algorithmus:
Bailey - Borwein - Plouffe - Formel
für Pi mit n = 5, n=1.000 (k=20), n=(10^9-1) und n=(10^10-1)
Diese Seite gibt es in zwei (drei)Versionen:
index.html (70k): HTML 4.0, "reines ASCII"
index.xhtml (210k): XHTML 1.1 plus (Presentation) MathML 2.0
index.xml: XHTML 1.1 plus (Content) MathML 2.0
# Worum geht's? | # Literatur | # Bezeichnungsweise | # Vorwort | # BBP | # Effizient modulo | # Vektor-, Parallel-isierung | # Fortran90 | # Weitere Links
Berechnung einer beliebigen (hier in diesen Beispielen die sechste, (teilweise) die 1.001ste, die einmilliardste und die zehnmilliardste) hexadezimalen Nachkommastelle von Pi ohne Kenntnis der vorhergehenden Ziffern.
Diese page ist aus der Absicht heraus entstanden, das Ganze selber genau verstehen zu wollen.
http://www.cecm.sfu.ca/personal/pborwein/
http://www.cecm.sfu.ca/personal/pborwein/PISTUFF/FORTRAN
http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf
Bezeichnungsweise:
(mathematisch) "hoch" (Exponentiation, manchmal auch mit ** symbolisiert)
"mal" (Mulitplikation)
bei einem Argument: Gauss-Klammer, Ganzzahlfunktion
bei mehreren Argumenten: Array, "Vektor" (Argumente sind Elemente/Komponenten)
x mod(ulo) y
Rest bei der Division von x durch y
bei nicht ganzen Zahlen: = x - [x/y]*y
Spezialfall y=1: x mod 1 = Nachkommateil von x = x - [x]
a ≡ b(m)
a "ist kongruent" b modulo m, d.h. a hat den gleichen Rest wie b bei Division durch m
innerhalb einer Zahl: Komma
innerhalb eines Arrays: Zahlen-Trennzeichen
Tausender- bzw. allg. Stellen-Trennzeichen (innerhalb einer Zahl)
Letztes (vor(her)iges) Ergebnis
Um z.B. die sechste (dezimale) Nachkommastelle von Pi (3,141.592.653.589.792...) zu berechnen, könnte folgendermaßen vorgegangen werden:
10^5 * Pi = 314.159,265.358.979.2...
mod(ulo) 1 = Nachkommateil: 0,265.358.979.2...
mal 10: 2,653.589.792...
Nachkommateil weglassen, [%]: 2
Das mag jetzt auf den ersten Blick künstlich aufwendig erscheinen, ist aber der grundlegende Gedanke für einen programmierbaren Algorithmus - bestehend u.a. aus wohldefinierten Rechenoperationen - , der von einem Computer abgearbeitet werden soll. Es könnte auch gleich mit 10^6 multipliziert werden, mod(ulo) 10 und dann Ganzzahl nehmen, aber gerade der mod 1 Teil ist für BBP wesentlich.
Für die sechste hexadezimale Nachkommastelle (Pi{16} = 3,243.F6A.888.5A...) sieht das dann so aus (im anschaulichen dezimalen System):
16^5 * Pi = 1.048.576 * 3,141.592.653.589.792... = 3.294.198,658.330.57...
mod(ulo) 1 = Nachkommateil: 0,658.330.57...
mal 16: 10,533.289.1...
[%]: 10 = A
Um solcherart z.B. die einmilliardste Stelle auszurechnen (was mit BBP schon gemacht wurde, s.a. Fortran90), müssen folgende zwei Probleme berücksichtigt werden:
Es ergäben sich riesige Zahlen bei Berechnung von 16^999.999.999 (* Pi), die, wenn überhaupt bewältigbar (u.a. Speicherproblem), zu extremen Rechenzeiten führten.
Damit das Ergebnis stimmt, müßte Pi schon am Anfang mit einer Milliarde Stellen genau bekannt sein , d.h. die gesuchte Ziffer (und auch alle vorhergehenden) müssen schon irgendwie anders ausgerechnet worden sein.
Problem I kann unter Ausnützung der Kommutativität der modulo-Rechnung mit den vier Grundrechnungsarten (und darausfolgend auch mit der Exponentiation) gelöst werden. Und zwar wird statt mit 16^5 nur mit 16 multipliziert und dafür Punkt 2 und 3 des obgenannten Algorithmus 6mal wiederholt (% bedeutet das Ergebnis der vorigen Rechnung/Zeile):
Pi mod 1	= 0,141.592.653.589
% * 16 = 2,265.482.457.436
% mod 1	= 0,265.482.457.436
% * 16 = 4,247.719.318.987
% mod 1	= 0,247.719.318.987
% * 16 = 3,963.509.103.793
% mod 1	= 0,963.509.103.793
% * 16 = 15,416.145.660.689 (15=F)
% mod 1	= 0,416.145.660.689
% * 16 = 6,658.330.571.034
% mod 1	= 0,658.330.571.034
% * 16 = 10,533.289.136.557
Beachte, daß die Zahlen vor dem Komma nach der Multiplikation mit 16 (Punkt 3) jeweils die hexadezimalen Ziffern darstellen: Pi{hex} = 3,243F6A...! (Wenn dieses Schema im dezimalen System - immer 10 statt 16 nehmen - ausgeführt wird, ist dies anschaulich klar.)
Der Vorteil dabei ist, daß immer nur Zahlen kleiner als 16 auftreten! Die zwei Nachteile: Erstens braucht mensch jetzt erst wieder alle vorherigen Ziffern (bzw. braucht nicht, aber erhält sie als Zwischenergebnisse, was irgendwie aufs gleiche herauskommt) und zweitens dauert die Rechnung bei z.B einer Milliarde Wiederholungen auch nicht gerade kurz. Zweiteres werden wir weiter unten behandeln, nachdem wir uns Problem II angeschaut haben:
Hier ist der Punkt, wo der BBP-Algorithmus ins Spiel kommt:
BBP (Bailey - Borwein - Plouffe):
Anstatt Pi wird ein numerisch äquivalenter Ausdruck gemäß folgender Formel verwendet (s.a.bbp1.png):
---   1    (   4       2       1       1   )
Pi = >   ---- * ( ----- - ----- - ----- - ----- )     :§1:
--- 16^k   ( 8*k+1   8*k+4   8*k+5   8*k+6 )
Mit der generellen Identität:
Z(ähler)	       Rest       Z mod N   Z mod N
Bruch = -------- = Ganzzahl + ---- = G + ------- ≡ ------- (1)
N(enner)	        N	     N	       N
führt dies zu (vgl. oben: um die 6. (n+1) Ziffer zu erhalten, wurde mit 16 ^ 5 (n) multipliziert) (s.a.bbp2.png):
[  n                      oo           ]
[ --- 16^(n-k) mod 8k+1   --- 16^(n-k) ]
16^n * Pi ≡ 4 * [ >   ----------------- + >   -------- ] -     :§2:
[ ---       8k+1          ---   8k+1   ]
[ k=0                    k=n+1         ]
[ --- 16^(n-k) mod 8k+4   --- 16^(n-k) ]
- 2 * [ >   ----------------- + >   -------- ] -
[ ---       8k+4          ---   8k+4   ]
[ --- 16^(n-k) mod 8k+5   --- 16^(n-k) ]
- 1 * [ >   ----------------- + >   -------- ] -
[ ---       8k+5          ---   8k+5   ]
[ --- 16^(n-k) mod 8k+6   --- 16^(n-k) ]
- 1 * [ >   ----------------- + >   -------- ] (1)
[ ---       8k+6          ---   8k+6   ]
Beachte, daß die Aufspaltung in jeweils zwei Summen (k=0..n und k=n+1..oo) nicht die Grenze zwischen Vor- und Nach-kommateil der Zahl beschreibt, sondern zwei Bereiche trennt, die verschieden berechnet werden können. Die Summen rechts sind wohl kleiner als eins, aber die linken tragen ebenso zum Nachkommateil bei.
Für n = 5 (gesucht wird die sechste Ziffer) ergeben sich die jeweils linken Summen gemäß folgender Tabelle:
Tabelle I: Teilsummen
16^(n-k)
16^(n-k) mod 8k+1
16^(n-k) mod 8k+4
16^(n-k) mod 8k+5
16^(n-k) mod 8k+6
0 16^5 1 0 4 0 5 1 6 4
1 16^4 9 7 12 4 13 3 14 2
2 16^3 17 16 20 16 21 1 22 4
3 16^2 25 6 28 4 29 24 30 16
4 16^1 33 16 36 16 37 16 38 16
5 16^0 41 1 44 1 45 1 46 1
s1, s2, s3, s4 0/1 + 7/9 + 16/17 + 6/25 + 16/33 + 1/41 =
2,468.192.977.116 0/4 + 4/12 + 16/20 + 4/28 + 16/36 + 1/44 =
1,743.362.193.362 1/5 + 3/13 + 1/21 + 24/29 + 16/37 + 1/45 =
1,760.629.139.939 4/6 + 2/14 + 4/22 + 16/30 + 16/38 + 1/46 =
1,967.467.086.689
Die dafür nötigen modulo-Rechnungen (Restklassenberechnungen) wurden laut Tabelle II durchgeführt, wobei wieder die Kommutativität ausgenützt wurde. Die fetten Zahlen sind jeweils die gesuchten Zwischenergebnisse für die obige Tabelle, die Rechnungen rechts daneben sind eigentlich nicht mehr nötig (daher kursiv gesetzt) aber der Anschaulichkeit halber auch angeführt. Die Rechnungen für 16^3 und 16^5 sind nur dort nicht kursiv, wo sie als Ergebnis für oben gebraucht werden, als Zwischenergebnis für andere 16^m Rechnungen sind sie nie nötig (siehe wieder weiter unten):
Tabelle II: 16^m modulo 8k+j
16^0 = 1 16^1 = 16 16^2 = 256 16^3 = 4.096 16^4 = 65.536 16^5 = 1.048.576
0 8k+1 1 1 ≡ 0 0 0 0 0 0
8k+4 4 1 16 ≡ 0 0 * 16 =0
0 ^ 2 = 0 0 * 16 = 0
0 ^ 3 = 0 0 * 16 = 0
0 ^ 4 = 0 0 * 16 = 0
0 ^ 5 = 0
8k+5 5 1 16 ≡ 1 1 * 16 = 16 ≡ 1
1 ^ 2 = 1 1 * 16 = 16 ≡ 1
1 ^ 3 = 1 1 * 16 = 16 ≡ 1
1 ^ 4 = 1 1 * 16 = 16 ≡ 1
1 ^ 5 = 1
8k+6 6 1 16 ≡ 4 4 * 16 = 64 ≡ 4
4 ^ 2 = 16 ≡ 4 4 * 16 = 64 ≡ 4
4 * 4 = 16 ≡ 4
4 ^ 3 = 64 ≡ 4 4 * 16 = 64 ≡ 4
4 ^ 2 = 16 ≡ 4
4 ^ 4 = 256 ≡ 4 4 * 16 = 64 ≡ 4
4 ^ 5 = 1.024 ≡ 4
1 8k+1 9 1 16 ≡ 7 7 * 16 = 112 ≡ 4
7 ^ 2 = 49 ≡ 4 4 * 16 = 64 ≡ 1
4 * 7 = 28 ≡ 1
7 ^ 3 = 343 ≡ 1 1 * 16 = 16 ≡ 7
4 ^ 2 = 16 ≡ 7
7 ^ 4 = 2.401 ≡ 7 7 * 16 = 112 ≡ 4
7 * 7 = 49 ≡ 4
7 ^ 5 = 16.807 ≡ 4
8k+4 12 1 ≡ 4 4 * 16 = 64 ≡ 4
8k+5 13 1 ≡ 3 3 * 16 = 48 ≡ 9
3 ^ 2 = 9 9 * 16 = 144 ≡ 1
9 * 3 = 27 ≡ 1
3 ^ 3 = 27 ≡ 1 1 * 16 = 16 ≡ 3
9 ^ 2 = 81 ≡ 3
3 ^ 4 = 81 ≡ 3 3 * 16 = 48 ≡ 9
3 ^ 5 = 243 ≡ 9
8k+6 14 1 ≡ 2 2 * 16 = 32 ≡ 4
2 ^ 2 = 4 4 * 16 = 64 ≡ 8
2 ^ 3 = 8 8 * 16 = 128 ≡ 2
8 * 2 = 16 ≡ 2
4 ^ 2 = 16 ≡ 2
2 ^ 4 = 16 ≡ 2 2 * 16 = 32 ≡ 4
8 * 4 = 32 ≡ 4
2 ^ 5 = 32 ≡ 4
2 8k+1 17 1 16 ≡ 1 1 * 16 = 16 16 * 16 = 256 ≡ 1
1 ^ 2 = 1 1 * 16 = 16
8k+4 20 1 16 ≡ 16 16 * 16 = 256 ≡ 16 16 * 16 = 256 ≡ 16
16 ^ 2 = 256 ≡ 16 16 * 16 = 256 ≡ 16
8k+5 21 1 16 ≡ 4 4 * 16 = 64 ≡ 1 1 * 16 = 16
4 ^ 2 = 16 16 * 16 = 256 ≡ 4
8k+6 22 1 16 ≡ 14 14 * 16 = 224 ≡ 4 4 * 16 = 64 ≡ 20
14 ^ 2 = 196 ≡ 20 20 * 16 = 320 ≡ 12
4 * 14 = 56 ≡ 12
3 8k+1 25 1 16 ≡ 6 6 * 16 = 96 ≡ 21 21 * 16 = 336 ≡ 11
6 ^ 2 = 36 ≡ 11 11 * 16 = 176 ≡ 1
21 * 6 = 126 ≡ 1
8k+4 28 1 16 ≡ 4 4 * 16 = 64 ≡ 8 8 * 16 = 128 ≡ 16
8k+5 29 1 16 ≡ 24 24 * 16 = 384 ≡ 7 7 * 16 = 112 ≡ 25
24 ^ 2 = 576 ≡ 25 25 * 16 = 400 ≡ 23
7 * 24 = 168 ≡ 23
8k+6 30 1 16 ≡ 16 16 * 16 = 256 ≡ 16 16 * 16 = 256 ≡ 16
4 8k+1 33 1 16 ≡ 25 25 * 16 = 400 ≡ 4 4 * 16 = 64 ≡ 31
25 ^ 2 = 625 ≡ 31 31 * 16 = 496 ≡ 1
4 * 25 = 100 ≡ 1
8k+4 36 1 16 ≡ 4 4 * 16 = 64 ≡ 28 28 * 16 = 448 ≡ 16
28 * 4 = 112 ≡ 4
8k+5 37 1 16 ≡ 34 34 * 16 = 544 ≡ 26 26 * 16 = 416 ≡ 9
34 ^ 2 = 1156 ≡ 9 9 * 16 = 144 ≡ 33
26 * 34 = 884 ≡ 33
8k+6 38 1 16 ≡ 28 28 * 16 = 448 ≡ 30 30 * 16 = 480 ≡ 24
28 ^ 2 = 784 ≡ 24 24 * 16 = 384 ≡ 4
30 * 28 = 840 ≡ 4
5 8k+1 41 1 16 ≡ 10 10 * 16 = 160 ≡ 37 37 * 16 = 592 ≡ 18
10 ^ 2 = 100 ≡ 18 18 * 16 = 288 ≡ 1
37 * 10 = 370 ≡ 1
8k+4 44 1 16 ≡ 36 36 * 16 = 576 ≡ 4 4 * 16 = 64 ≡ 20
36 ^ 2 = 1.296 ≡ 20 20 * 16 = 320 ≡ 12
4 * 36 = 144 ≡ 12
8k+5 45 1 16 ≡ 31 31 * 16 = 496 ≡ 1 1 * 16 = 16
31 ^ 2 = 961 ≡ 16 16 * 16 = 256 ≡ 31
1 * 31 = 31
8k+6 46 1 16 ≡ 26 26 * 16 = 416 ≡ 2 2 * 16 = 32
26 ^ 2 = 676 ≡ 32 32 * 16 = 512 ≡ 6
2 * 26 = 52 ≡ 6
16^0 = 1 16^1 = 16^0 * 16 = 16 16^2 = 16^1 * 16 = (16^1)^2 = 256 16^3 = 16^2 * 16 = 16^2 * 16^1 = (16^1)^3 = 4.096 16^4 = 16^3 * 16 = 16^3 * 16^1 = (16^2)^2 = (16^1)^4 = 65.536 16^5 = 16^4 * 16 = 16^4 * 16^1 = 16^3 * 16^2 = (16^1)^5 = 1.048.576
Wie anhand dieser Tabelle (und speziell deren Fußzeile) ersichtlich, gibt es verschiedene Möglichkeiten, um ein modulo-Ergebnis für 16^m auszurechen, ohne mit der (großen) Zahl 16^m selbst rechnen zu müssen (% ist wieder das Ergebnis der vorherigen/obigen Rechnung/Zeile, %% noch eine Rechnung/Zeile davor/darüber):
Ergebnis vorige Spalte "16^(m-1)"
% * 16
% mod 8k+j
Ergebnis voriger Spalte "16^(m-1)"
Ergebnis von Spalte "16^1"
%% * %
Verallgemeinerung von II: Ist m=a+b, also 16^m = 16^(a+b) = 16^a * 16^b, dann:
Ergebnis Spalte "16^a"
Ergebnis Spalte "16^b"
Es kann hier verschiedene Möglichkeiten für die Zerlegung geben, z.B.: m = 5 = 4 + 1 = 3 + 2.
Ist m=a*b, also 16^m = 16^(a*b) = (16^a)^b, also eine (ganzzahlige (b)) Potenz von 16^a, dann:
% ^ b
Es kann hier verschiedene Möglichkeiten für die Zerlegung geben, z.B.: m = 12 = 6 * 2 = 4 * 3 = 3 * 4 = 2 * 6.
Für's Kopfrechen ist die beste (einfachste) Variante jene, wo die kleinsten Zahlen auftreten. Für den Computer (bzw. die Rechenzeit) wäre aber die Abfrage aller Möglichkeiten, um zu sehen, wo die kleinsten Zahlen auftreten, auch nicht optimal. Es sei hier wiederum auf weiter unten verwiesen für die im Allgemeinen effizienteste Art der Berechnung.
Auf jeden Fall treten immer nur Zahlen kleiner (8k+j)^2 auf! (Bzw. kann sogar auf ((8k+j)/2)^2 = ((8k+j)^2)/4, also auf ein Viertel reduziert werden, wenn man negative Zahlen dazunimmt; z.B.: 48 ≡ 9 ≡ -4 (13).)
Das Ergebnis für die linken Summen von §2 ist nun laut Tabelle I :
4 * 2,468.192.977.116... -
- 2 * 1,743.362.193.362... -
- 1 * 1,760.629.139.939... -
- 1 * 1,967.467.086.689... =
=     2,657.951.295.114...    :§3:
macht modulo 1:
≡ 0,657.951.295.114... (1)     :§4:
(vgl. den Wert mit der 6. Wiederholung von Punkt 2 oben!)
noch ein Mal mit 16 multipliziert:
% * 16 = 10,527.220.721.837...     :§5:
Der Vorkommateil ist nun die 6. hexadezimale Ziffer von Pi:
[%] = 10{16} = A :§6:
Versuchen wir auch noch die nächsten Ziffern zu errechnen:
§5 mod 1 = 0,527.220.721.837...
% *  16 = 8,435.531.549.394...
% mod 1 = 0,435.531.549.394...
% *  16 = 6,968.504.790.309...
D.h. die nächste Ziffer 8 kommt noch richtig heraus, die übernächste ist schon falsch. Korrigierien wir §3 um die jeweils ersten Terme der rechten Summen in §2, k = n+1 = 6: (s.a.bbp3.png)
4              2              1              1
------------ - ------------ - ------------ - ------------ =     :§7:
16 * (8*6+1)   16 * (8*6+4)   16 * (8*6+5)   16 * (8*6+6)
= 0,000.361.541.972...
Damit wird §4 zu (vgl. wieder 6. Wh. oben):
0,658.312.837.086...
und die nächsten Ziffern:
% * 16 = 10,533.005.393.384...
% mod 1 = 0,533.005.393.384...
% *  16 = 8,528.086.294.157...
% mod 1 = 0,528.086.294.157...
% *  16 = 8,449.380.706.525...
% mod 1 = 0,449.380.706.525...
% *  16 = 7,190.091.304.409...
also eine richtige Ziffer mehr.
Der nächste Korrektur-Term laut §2 wäre dann: k = n+2 = 7 (s.a.bbp4.png):
4                2                1                1
-------------- - -------------- - -------------- - -------------- =
16^2 * (8*7+1)   16^2 * (8*7+4)   16^2 * (8*7+5)   16^2 * (8*7+6)
= 0,000.016.873.556...
(das ist ca. ein 21-tel des ersten Korrekturterms,) womit mindestens eine weitere richtige Ziffer extrahierbar wird. D.h. durch entsprechend viele Terme der rechten Summe in §2 erhält mensch "beliebig" viele richtige Ziffern (hexadezimale Nachkommastellen) ab der gewünschten, ohne die vorherigen zu kennen (bzw. kennen zu müssen); vorausgesetzt natürlich, es wurde dauernd mit der entsprechenden Rechengenauigkeit (Stellenanzahl) gerechnet!
Nun aber zum oftmals angekündigten Thema:
Effizientes modulo Rechnen:
Beispiel für n = 1.000, k = 20 und s1 (8k+1):
Gesucht wird also der Ausdruck
16^(n-k) mod 8k+1 = 16 ^ 980 mod 161
16^980 ist ungefähr 10^1.180 also eine Zahl mit über 1.180 Ziffern, welche zwar durchaus noch handhabbar wäre für Computer, aber bei n = 1 Milliarde wären wir weit jenseits. Ich nehme hier 1.000 statt 1.000.000.000, damit der Umfang der page im Rahmen bleibt; es geht ja nur darum, das Prinzip zu zeigen:
Wir könnten nun also - beginnend bei 1 - 980 mal mit 16 multiplizieren und jeweils modulo 161 rechnen, das wären 980 Multiplikationen und ebensoviele modulo Rechnungen. Mit der hier nun vorzustellenden binären Methode reduziert sich der Aufwand auf jeweils maximal 10 (plus eine Binärzerlegung):
-1.: 16^  0 = . . . . . . . . . . . . . . 1 ≡   1 (161)
0.: 16^  1 =  . . . . . . . . . . . . . 16 ≡  16 (161)
1.: 16^  2 = (16^  1) ^ 2 ≡  16^2 =    256 ≡  95 (161)
2.: 16^  4 = (16^  2) ^ 2 ≡  95^2 =  9.025 ≡   9 (161)
3.: 16^  8 = (16^  4) ^ 2 ≡   9^2 =     81 ≡  81 (161)
4.: 16^ 16 = (16^  8) ^ 2 ≡  81^2 =  6.561 ≡ 121 (161)
5.: 16^ 32 = (16^ 16) ^ 2 ≡ 121^2 = 14.641 ≡ 151 (161)
6.: 16^ 64 = (16^ 32) ^ 2 ≡ 151^2 = 22.801 ≡ 100 (161)
7.: 16^128 = (16^ 64) ^ 2 ≡ 100^2 = 10.000 ≡  18 (161)
8.: 16^256 = (16^128) ^ 2 ≡  18^2 =    324 ≡   2 (161)
9.: 16^512 = (16^256) ^ 2 ≡   2^2 =      4 ≡   4 (161)
D.h. also wiederholtes (hier 10 mal) Quadrieren und modulo 161 Rechnen.
Die nächste Zeile (10.) wäre 16^1.024, das ist aber schon mehr als die benötigten 16^980, daher können wir hier aufhören. Alle beteiligten Zahlen sind kleiner als 25.921 = 161^2. (Wenn im Vorhinein bekannt ist, welches maximale n verwendet wird, können die Variablen entsprechend speicherplatz- und in Folge auch rechenzeit-schonend deklariert werden.)
Als nächstes brauchen wir die binäre Darstellung von 980:
980{2} = 11.1101.0100
|| || |  |
980 = 2^9 + 2^8 + 2^7 + 2^6 + 2^4 + 2^2 =
= 512 + 256 + 128 +  64 +  16 +   4
Folglich (unter Benutzung der vorigen Ergebnisse):
16^980 = 16 ^ (512 +    256 +    128 +    64 +    16 +    4) =
=    16^512 * 16^256 * 16^128 * 16^64 * 16^16 * 16^4  ≡
≡         4 *      2 *     18 *   100 *   121 *    9 (161)
Die letzte Rechnung wird (v.a. bei größeren Zahlen) zweckmäßigerweise auch wieder sukzessive durchgeführt (damit bleiben die Zwischenergebnisse wieder unter 161^2):
4 *   2 =      8 ≡   8 (161)
8 *  18 =    144 ≡ 144 (161)
144 * 100 = 14.400 ≡  71 (161)
71 * 121 =  8.591 ≡  58 (161)
58 *   9 =    522 ≡  39 (161)
Auch hier wären maximal 10 Zyklen - je nachdem, wieviele 1er in der binären Darstellung von n-k besetzt sind - zu rechnen, wenn n=1.000, weil: log(1.000)/log(2) = 9,965.784...
16^980 ≡ 39 (161)
Das ging ja recht flott und war eigentlich gar nicht so schwierig, oder?
Dieses eine Beispiel war für k=20 (das Ergebnis ist: 39/161 = 0,242.236.025...), insgesamt (um wirklich die 1.001ste Hex-Ziffer von Pi zu berechnen) muß das gleiche für alle k von 0 bis n (=1.000; d.h. mensch muß maximal modulo 8k+1 = 8.001 rechnen) gemacht werden; und dann das Ganze jeweils auch für s2 (8k+4), s3 (8k+5) und s4 (8k+6); die vier Teilergebnisse richtig summieren, mod 1, mal 16, und der Ganzzahlteil davon ist das Ergebnis, die gesuchte Hex-Ziffer. Das ist durchaus kein kleiner Aufwand, hat mensch den aber einmal gemacht, dann ist die zusätzliche Arbeit für weitere Ziffern vergleichbar winzig, nur vier Ausdrücke müssen ausgewertet werden, nämlich die jeweils ersten Terme der rechten Summen in §2, um das Ergebnis der "Monster"-Rechnung zu verbessern (vgl. §7, k=n+1=1.001): (s.a.bbp5.png)
4                  2                  1                  1
---------------- - ---------------- - ---------------- - ---------------- =
16 * (8*1.001+1)   16 * (8*1.001+4)   16 * (8*1.001+5)   16 * (8*1.001+6)
= 0,000.000.014.608.4... = ~ 16^(-6,5)
woraus auch ersichtlich wird, daß diese Korrektur erst 6 Ziffern weiter zum Tragen kommt, also ohne die Korrektur schon fünf bis sechs Ziffern richtig sind.
Bei n = 1 Milliarde sind mit der binären Methode jeweils maximal 30 (statt 1 Mia.) Zyklen zu rechnen (log(10^9)/log(2) = 9*log(10)/log(2) = 29,897.352...), also eine ansehnliche Ersparnis.
Für n = 5 ist die Einsparung noch nicht so überwältigend, aber Tabelle II wird zu:
Tabelle III: 16^(2^i) modulo 8k+j (effizient) für n=5
n-k ({2})
0 8k+1 1 16 ≡ 0 0 0 0 * 0 = 0 0 / 1 = 0 5
8k+4 4 16 ≡ 0 0 ^ 2 = 0 0 ^ 2 = 0 0 * 0 = 0 0 / 4 = 0
8k+5 5 16 ≡ 1 1 ^ 2 = 1 1 ^ 2 = 1 1 * 1 = 1 1 / 5 = 0,2
8k+6 6 16 ≡ 4 4 ^ 2 = 16 ≡ 4 4 ^ 2 = 16 ≡ 4 4 * 4 = 16 ≡ 4 4 / 6 = 0,666.666.666...
1 8k+1 9 16 ≡ 7 7 ^ 2 = 49 ≡ 4 4 ^ 2 = 16 ≡ 7 7 7 / 9 = 0,777.777.777... 4
8k+4 12 16 ≡ 4 4 ^ 2 = 16 ≡ 4 4 ^ 2 = 16 ≡ 4 4 4 / 12 = 0,333.333.333...
8k+5 13 16 ≡ 3 3 ^ 2 = 9 9 ^ 2 = 81 ≡ 3 3 3 / 13 = 0,230.769.231...
8k+6 14 16 ≡ 2 2 ^ 2 = 4 4 ^ 2 = 16 ≡ 2 2 2 / 14 = 0,142.857.143...
2 8k+1 17 16 16 ^ 2 = 256 ≡ 1 16 * 1 = 16 16 / 17 = 0,941.176.471... 3
8k+4 20 16 16 ^ 2 = 256 ≡ 16 16 * 16 = 256 ≡ 16 16 / 20 = 0,8
8k+5 21 16 16 ^ 2 = 256 ≡ 4 16 * 4 = 64 ≡ 1 1 / 21 = 0,047.619.047...
8k+6 22 16 16 ^ 2 = 256 ≡ 14 16 * 14 = 224 ≡ 4 4 / 22 = 0,181.818.181...
3 8k+1 25 16 16 ^ 2 = 256 ≡ 6 6 6 / 25 = 0,24 2
8k+4 28 16 16 ^ 2 = 256 ≡ 4 4 4 / 28 = 0,142.857.143...
8k+5 29 16 16 ^ 2 = 256 ≡ 24 24 24 / 29 = 0,827.586.207...
8k+6 30 16 16 ^ 2 = 256 ≡ 16 16 16 / 30 = 0,533.333.333...
4 8k+1 33 16 16 16 / 33 = 0,484.848.484... 1
8k+4 36 16 16 16 / 36 = 0,444.444.444...
8k+5 37 16 16 16 / 37 = 0,432.432.432...
8k+6 38 16 16 16 / 38 = 0,421.052.632...
5 8k+1 41 1 1 / 41 = 0,024.390.243... 0
8k+4 44 1 1 / 44 = 0,022.727.272...
8k+5 45 1 1 / 45 = 0,022.222.222...
8k+6 46 1 1 / 46 = 0,021.739.130...
16^(5-k) mod 8k+j
(<- %) / 8k+j
Vektor-, Parallel-isierung:
Alle auszuführenden Rechenoperationen sind für j=1,4,5,6 (s1,s2,s3,s4) gleichartig, nur mit unterschiedlichen Zahlen, daher können sie auf Rechnern bzw. mit Software, welche das beherrschen/beherrscht (quasi-)parallel durchgeführt werden, indem statt mit einfachen Zahlen mit 4er-Vektoren bzw. 4er-Arrays gerechnet wird, wobei die mathematischen Operationen elementweise durchzuführen sind und nicht als Vektoroperationen, z.B.:
[1,2,3,4] + [10,20,30,40] = [1+10,2+20,3+30,4+40] = [11,22,33,44]
[4,3,2,1] * [10,20,30,40] = [4*10,3*20,2*30,1*40] = [40,60,60,40]
[1,2,3,4] ^ [ 4, 3, 2, 1] = [1^4 ,2^3 ,3^2 ,4^1 ] = [ 1, 8, 9, 4]
[1,2,3,4] ^ 2 = [1^2,2^2,3^2,4^2] = [ 1, 4, 9,16] ≡ [1,4,2,0] mod [5,6,7,8]
Mit dieser Notation erhalten wir nun schließlich die Tabelle für n = 1.000 und k = 20 (sowie ein paar weitere), wobei die laufende Berechnung der (Zwischen-)Ergebnisse auch gleich integriert ist, zu:
Tabelle IV: 16^(2^i) modulo 8k+j für n=1.000 (, k=20)
8k+[1,4,5,6]
n-k {2}
16^64
16^128
16^256
16^512
0 [1,4,5,6] [0,4,5,6] [0,0,1,4] [0,0,1,4] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] / [1,4,5,6] = [,,,]
1.000 {2} = 0 0 0 1 0 1 1 1 1 1
[,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,]
1 [9,12,13,14] [7,4,3,2] [4,4,9,4] [7,4,3,2] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] / [9,12,13,14] = [,,,]
999 {2} = 1 1 1 0 0 1 1 1 1 1
[7,4,3,2] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,]
2 [17,20,21,22] [16,16,16,16] [1,16,4,14] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] / [17,20,21,22] = [,,,]
998 {2} = 0 1 1 0 0 1 1 1 1 1
[1,16,4,14] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,]
3 [25,28,29,30] [16,16,16,16] [6,4,24,16] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] / [25,28,29,30] = [,,,]
997 {2} = 1 0 1 0 0 1 1 1 1 1
[16,16,16,16] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,]
4 [33,36,37,38] [16,16,16,16] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] / [33,36,37,38] = [,,,]
996 {2} = 0 0 1 0 0 1 1 1 1 1
[,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,]
5 [41,44,45,46] [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] ^ 2 = [,,,] ≡ [,,,] [,,,] / [41,44,45,46] = [,,,]
995 {2} = 1 1 0 0 0 1 1 1 1 1
[,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,] [,,,] * [,,,] = [,,,] ≡ [,,,]
... {2} =
20 [161,164,165,166] [16,16,16,16] [16,16,16,16] ^ 2 = [256,256,256,256] ≡ [95,92,91,90] [95,92,91,90] ^ 2 = [9.025,8.464,8.281,8.100] ≡ [9,100,31,132] [9,100,31,132] ^ 2 = [81,10.000,961,17.424] ≡ [81,160,136,160] [81,160,136,160] ^ 2 = [6.561,25.600,18.496,25.600] ≡ [121,16,16,36] [121,16,16,36] ^ 2 = [14.641,256,256,1.296] ≡ [151,92,91,134] [151,92,91,134] ^ 2 = [22.801,8.464,8.281,17.956] ≡ [100,100,31,28] [100,100,31,28] ^ 2 = [10.000,10.000,961,784] ≡ [18,160,136,120] [18,160,136,120] ^ 2 = [324,25.600,18.496,14.400] ≡ [2,16,16,124] [2,16,16,124] ^ 2 = [4,256,256,15.376] ≡ [4,92,91,104] [39,124,1,144] / [161,164,165,166] = [0,242.236.025..., 0,756.097.561..., 0,006.060.606..., 0,867.469.880...]
980 {2} = 0 0 1 0 1 0 1 1 1 1
[9,100,31,132] [9,100,31,132] * [121,16,16,36] = [1.089,1.600,496,4.742] ≡ [123,124,1,104] [123,124,1,104] * [100,100,31,28] = [12.300,12.400,31,2.912] ≡ [64,100,31,90] [64,100,31,90] * [18,160,136,120] = [1.152,16.000,4.216,10.800] ≡ [25,92,91,10] [25,92,91,10] * [2,16,16,124] = [50,1.472,1.456,1.240] ≡ [50,160,136,78] [50,160,136,78] * [4,92,91,104] = [200,14.720,12.376,8.112] ≡ [39,124,1,144]
... 512 {2} =
489 ... 744
511 ... 256 {2} =
745 ... 872
255 ... 128 {2} =
873 ... 936
127 ... 64 {2} =
937 ... 968
63 ... 32 {2} =
969 ... 984
31 ... 16 {2} =
985 ... 992
15 ... 8 {2} =
993 ... 996
7 ... 4 {2} =
997 ... 998
3 ... 2 {2} =
1 {2} =
0 {2} =
16^(2^0)
16^(2^1)
16^(2^2)
16^(2^3)
16^(2^4)
16^(2^5)
16^(2^6)
16^(2^7)
16^(2^8)
16^(2^9)
16^(n-k) mod 8k+j / 8k+j
(Die vollständige Tabelle ist ca. 1,7M groß!)
Die hellgraue Hinterlegung dient nur der Übersichtlichkeit.
Die Binärdarstellung von n-k läuft rückwärts, passend zu 16^(2^i) !
Die Zeilen für k=0...5 wurden deshalb tlw. mit Inhalt übernommen, um zu zeigen, daß die (16^(2^i) mod 8k+j)-Ergebnisse von Tab. III übernommen werden können. D.h. möchte man später für größere n rechnen, können die schon gerechneten Werte weiterverwendet werden, die ("dreieckige") Tabelle muß "nur" nach rechts (und nach) unten erweitert werden.
Mit steigendem k muß nicht mehr zu so hohen Potenzen von 16 gerechnet werden (n-k), dafür werden die beteiligten Zahlen (wegen der steigenden 8k+j) größer.
Auf die öfters vorkommenden (zyklischen) Wiederholungen in den modulo Rechnungen (z.B. auch bei k=20) wird hier nicht näher eingegangen.
Programmiertechnisch werden durch die laufende Mitberechnung nur drei Vierer-Arrays ("VA") als Variablen gebraucht (abgesehen natürlich von den Zählschleifen-Variablen und anderem Drumherum):
VA1: enthält die (16^(2^i) mod 8k+j)-Werte; von einer Spalte zur nächsten (i-Schleife) lautet die Berechnung:
VA1 = VA1^2 mod (8*k+[1,4,5,6])
VA2: enthält die Zwischenergebnisse gemäß 3. Subzeile innerhalb einer k-Zeile in Tab. IV:
if i-te Binärstelle von n-k==1 then:
VA2 = VA1*VA2 mod (8*k+[1,4,5,6])
VA3: enthält die laufend aufsummierten Zwischenergebnisse gemäß der letzten Spalte in Tab. IV und wird nicht innerhalb der i-Schleife, sondern an derem Ende, in der k-Schleife neu berechnet:
VA3 = ( VA3 + VA2 / {8*k+[1,4,5,6]}  ) mod 1
mod 1 bewirkt, daß nur der nötige Nachkommateil mitgeschleppt wird und somit die Rechengenauigkeit in den letzten Stellen nicht (so stark) durch Aufsummierungseffekte verloren geht. Achtung: so programmieren, daß (8k+j) / VA2 nicht als Integer-Division ausgeführt wird!
Das Endergebnis lautet dann:
Gesuchte Ziffer = [ 16 * ( [4,-2,-1,-1]*VA3 mod 1 ) ]
bbp.f90.txt:
Ein ganz kurzes (1.5kB) Fortran90 Programm zum Demonstrieren, wie kompakt das sein kann. Es funktioniert bis n=5.791, weil: k geht von 0 bis n, 8k+j wird daher maximal 8n+6, der entsprechende modulo-Rest maximal 8n+5. Da innerhalb der Berechnung nur (Integer-)Zahlen kleinergleich (8k+j)^2 auftreten, muß für die größte, vom Programm/Compiler verwendbare Integerzahl I gelten: I >= (8n+5)^2, oder: n <= (sqrt(I)-5)/8. sqrt() bedeutet Wurzel(), die größte Zahl wird in Fortran90 mit huge() abgefragt.
bbpI2.f90.txt:
Das Ganze mit doppelt langen Integer - da geht's dann gleich bis n=379.625.061 - und ein paar zusätzlichen Features. Achtung: die "final sum" ist nicht so genau, wie durch den 1. Korrekturterm angegeben wird und das maximale n (für richtige Ergebnisse) ist auch kleiner!: Die Summe (VA3, arr3) hat eine Genauigkeit von 7 Ziffern, d.h. bei jeder Summation ist ein Fehler von bis zu 5*10^-8 drin, macht insgesamt (n+1)*5*10^-8. Z.B. n=99.999: Final sum = 0.3256860, 1st corr. = 1.464833E-12, Result: 100000 th hex digit of pi = 5. Obwohl der Korrekturterm von der Größe 10^-12 ist, sind von der final-sum nur die ersten 2 Ziffern signifikant! Und da sind wir auch schon bald beim n-Limit: damit die hex-Ziffer verläßlich ist, darf der Fehler nicht mehr als 8*16^-2 betragen, das ergibt: n < 8*10^8/(5*16^2) = 625.000. Um also die doppelt langen Integer wirklich nutzen zu können, müssen wir auch bei den Real was zulegen:
bbpI2R2.f90.txt:
Zusätzlich doppelt lange Real, beim Ergebnis wird auch der geschätzte Fehler (estimated error) angegeben.
bbpI2R4.f90.txt:
Vierfach lange Real, beim Ergebnis werden auch noch weitere Ziffern (Anzahl je nach Korrektionsterm) angegeben, wobei die letzte nicht mehr stimmen muß. Die Ziffern in Klammer sind unsicher aufgrund des geschätzten (Summierungs-)Fehlers.
bbpR4R4.f90.txt:
Ein zweites Mal vierfach lange Real, in Ermangelung von vierfach langen Integer. Natürlich dürfen diese Pseudo-Integer nicht bis zur Real-Grenze (~1.1897E+4932) ausgereizt werden, sondern als größte Zahl ist jene zu nehmen, die sich nicht mehr ändert, wenn 1 addiert wird, aber schon ändert, wenn 1 subtrahiert wird (~1.03846E+0034) Als maximales n erhalten wir hier 12.738.103.345.051.544. Es werden so viele Korrektionsterme wie nötig dazugerechnet um die volle Zahlengenauigkeit auszunutzen, allenfalls wieder durch den geschätzten Fehler begrenzt.
Für die einmilliardste (und folgenden) Stelle(n) erhalte ich damit (nach ca. 60 Stunden Rechenzeit, welche proportional zu n*log(n) ist):
| 'bbpR4R4.f90':							   |
|   Fortran90 program implementing the BBP algorithm			   |
|     for computing hex digits of Pi according to			   |
|     http://www.unet.univie.ac.at/~a8727063/Science/BBP/		   |
| Integer Kind = 8|16 , Real Kind = 16					   |
| max n = 12738103345051544 = 1.27E+16 (6.25E+29)			   |
exit program with n < 0
n = ? : 999999999
Final sum =   0.5216268016085995792268670859006806
Est. err. =    5.000000000000000000000000000000001E-0023
1st corr. =    1.464843748913574219459533606131920E-0020
Corr. sum =   0.5216268016085995792424920858870090	(   13 corr.s)
Result: 1000000000 th and following hex digits of pi =
8 5 8 9 5 5 8 5 A 0 4 2 8 B 5 6 ( 4 0 8 4 E 7 4 9 E 4 E 8 0 )
n = ? : -1
Für die zehn-milliardste (nach 2.327.073,607 Sekunden = 27 Tage Rechenzeit):
n = ? : 9999999999
Final sum =   0.5707466468480209374959415108034570
Est. err. =    5.000000000000000000000000000000001E-0022
1st corr. =    1.464843749891357421882095580551268E-0022
Corr. sum =   0.5707466468480209374960977608034433	(   11 corr.s)
Result: 10000000000 th and following hex digits of pi =
9 2 1 C 7 3 C 6 8 3 8 F B 2 B ( 6 2 2 3 6 3 0 F 3 8 5 9 A 0 )
http://www.lupi.ch/PiSites/pibinhexdec.html:
Ihr Browser kann keine iframes anzeigen, bitte benützen Sie den obigen link.
http://www.google.at/search?q=Bailey+Borwein+Plouffe
Sollte ich mich irgendwo verschrieben haben, oder noch etwas grob unklar sein, bitte die unten angegebene Adresse verwenden.
sa.13.dez.2k+8, © & @ (= i-mehl): Johannes Nendwich

References: §2

§5
 §3
 §2
 §4
 §2
 §2
 §2
 §7