Leuk dat dit topic weer boven komt borrelen. Ik heb toen ik mijn promotieonderzoek ff niet meer zag zitten, me een tijdje verdiept in hoe/wat van de praktische berekening van PI. Alles wat ik toen heb uitgezocht, heb ik in LaTeX vorm gegoten bij wijze van "oefen" publicatie.
Ik geef het hierbij aan GoT met de hoop dat iemand dit nog even door LaTeX kan halen en in een wat gemakkelijker formaat ergens beschikbaar kan stellen voor alle pi-tweakers (in Word, HTML?).
Thanx,
Martin
\documentstyle[11pt]{article}
\title{ {\huge De Berekening van $\pi$ in de Praktijk} }
\author{ \\ \\
{Dr. Martin E. Los} \\ \\ \\ \\
}
%\begin{document}
\maketitle
\begin{abstract}
{\noindent
In dit artikel laat ik zien, hoe je in de praktijk een
groot aantal decimalen van $\pi$ kan berekenen.
Allereerst geef ik een uitgebreid historisch overzicht
van de jacht op decimalen van $\pi$ door de eeuwen heen.
Vervolgens worden de snelle Borwein algorithmes
besproken, die gebruikt worden in alle moderne berekeningen.
Tenslotte geef ik praktische aanwijzingen voor het oplossen
van de problemen die je tegen komt wanneer je deze Borwein
algorithmes wilt implementeren.
Hopelijk is dit artikel een compact naslagwerk over $\pi$.
Wie nog meer details wil weten, verwijs ik naar de omvangrijke
literatuurlijst aan het eind van dit artikel.
}
\end{abstract}
\begin{document}
\newpage
\section{Een Historisch Overzicht}
De eerste~\cite{BECKMANN} die $\pi$ interessant~\cite{SAGAN}
genoeg vond was Archimedes. Hij probeerde $\pi$ te benaderen door
de oppervlaktes te berekenen van gelijkzijdige polygons met een
groeiend aantal zijden. Triest genoeg kwam hij niet verder dan
$\frac{6336}{2017.25} < \pi < \frac{14688}{4673.5}$
( $3.1405 < \pi < 3.1428$ ).
Dit krijg je met 96 gelijkzijdige polygons (probeer dit niet thuis!).
Archimedes gooide er het bijltje bij neer en vond dat $\frac{22}{7}$
in de praktijk voldeed.
In begin 1600 berekende een 'idiot savant' genaamd Ludolph von Ceulen
35 decimalen van $\pi$ met behulp van Archimedes' methode.
Gelukkig vond Gregory in 1671 de arctangens serie (of Taylor expansie): \\
$\arctan x = x - \frac{x^{3}}{3} + \frac{x^{5}}{5} - \frac{x^{7}}{7} + \cdots $ ,
zodat $\frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7}
+ \cdots $ \\
\noindent
Deze reeks en variaties die ietwat sneller convergeren, leidden in de
komende 300 jaar tot een paar vermeldenswaardige berekeningen.
In 1706 deed John Machin met $\pi = 16 \arctan (\frac{1}{5}) - 4 \arctan (
\frac{1}{239})$ een succesvolle gooi naar de honderd decimalen. Johann Dase
berekende in 1844 met
$\pi = 4 \arctan (\frac{1}{2}) + 4 \arctan (\frac{1}{5}) +
4 \arctan (\frac{1}{8})$ tweehonderd decimalen en tenslotte produceerde
William Shanks met dezelfde reeks in 1873 zelfs 707 goede decimalen.
Metius, een Nederlands wiskundige rond 1700, kwam met de eerste
praktische bijdrage. Hij vond de (4-de orde kettingbreuk~\cite{CIJSOUW})
benadering voor $\pi$ die ook nog gemakkelijk te onthouden is.
Schrijf de oneven numbers 1, 3 en 5 elk tweemaal achter elkaar : 113355.
Snijd vervolgens dit getal in de twee delen 113 en 355, deel 355 door 113
en voil\`{a} het antwoord $3.14159292$ is bijna $3,141592654 \dots$
In 1949 werd de eerste computerpoging gewaagd door
Rietweisner~\cite{RIETWEISNER} op de beroemde kamervullende ENIAC en hij
vond 2037 decimalen gebruik makend van Machin's reeks.
Daniel Shanks (misschien wel familie van William, geen idee ...)
berekende in 1962~\cite{SHANKS} samen met John Wrench en hun beider vriend
de IBM 7090 (die het meeste werk deed), ruim $100.000$ decimalen gebruik
makend van
$\pi = 24 \arctan (\frac{1}{8}) + 8 \arctan (\frac{1}{57}) +
4 \arctan (\frac{1}{239})$.
De echte doorbraak kwam in 1976 toen Brent~\cite{BRENT} en
Salamin~\cite{SALAMIN} onafhankelijk van elkaar een algorithme vonden
om $\pi$ te berekenen dat een kwadratische convergentie heeft.
Het is vreemd dat het zolang duurde voordat dit algorithme gevonden werd,
want de noodzakelijke relaties waren allang bekend bij Gauss~\cite{GAUSS}
en Legendre. Tevens had in het begin van deze eeuw de Indiase wiskundige
Ramanujan~\cite{RAMANUJAN} in zijn beroemde maar onleesbare 'notebooks'
ook al een kwadratisch convergerende formule neergekalkt, te weten :
\\[
\frac{1}{\pi} = \frac{1}{4} \left( \frac{1123}{882} -
\frac{22583}{882^{3}} \frac{1}{2} \frac{1 \cdot 3}{4^{2}} +
\frac{44043}{882^{5}}
\frac{1 \cdot 3}{2 \cdot 4} \frac{1 \cdot 3 \cdot 5 \cdot 7}{4^{2} \cdot 8^{2}} - \cdots
\right)
\]
De eerste factor in de teller gaat als $1123+21460 k$, met
$k=0,1,2,3,\ldots$
In de begin jaren tachtig hebben de Borwein's~\cite{BORAGM,BORPI,NEWMAN}
(2 broertjes neem ik aan)
de methode van Brent en Salamin gegeneraliseerd, hetgeen leidde tot een
een klasse van algorithmes voor het berekenen van $\pi$. Deze algorithmes
hebben de eigenschap dat ze bij elke iteratie $n$ keer meer decimalen
van $\pi$ op te leveren. Hoewel $n$ heel groot kan zijn, worden in de
praktijk alleen de algorithmes met lage $n$'s gebruikt omdat de hogere
niet effectiever zijn in een computer implementatie wegens hun hogere
complexiteit.
Gebruik makend van deze algorithmes berekenden Kanada en Tamura in
1983~\cite{KANADA} al 16 miljoen decimalen en David Bailey~\cite{BAILEYPI}
in 1986 $29.360.000$ decimalen.
In 1987 rapporteerde Kanada dat hij meer dan 134 miljoen decimalen
had berekend en zover ik weet gaat men nog altijd door ... \\
\noindent
Nu vraag je je misschien af waarom 'serieuse' mensen dit enorme aantal
decimalen van $\pi$ nodig hebben?
Allereerst is het nog niet bewezen dat de decimalen van $\pi$ willekeurig
zijn, hoewel wel bewezen is dat $\pi$ irrationeel is (Lambert, 1766).
Gegenereerde decimalen kunnen dus, omdat ze tot nu toe alle willekeurigheids testen doorstaan hebben, gebruikt worden om grote en niet repeterende
willekeurige getallenrijen te genereren.
Verder vormt dit soort berekeningen een prima test voor de
betrouwbaarheid van een computer. Het gaat immers om triljoenen
arithmetische operaties die tot een specifieke rij getallen moet leiden.
In alle latere computerpogingen is $\pi$ ter controle dan ook twee
keer berekenend met verschillende algorithmes, omdat de kans op fouten niet
meer verwaarloosbaar meer is bij berekeningen van deze omvang.
\section{De Borwein Algorithmes}
Ik zal hieronder de Borwein algorithmes met $n=2$ en $n=4$ introduceren.
Merk op dat je multi-precisie tools nodig hebt om ze te kunnen gebruiken,
vanaf het begin moet namelijk gerekend worden met het gewenste aantal
eind decimalen.\\
\noindent
Voor het kwadratisch ($n=2$) convergerende algorithme voor $\pi$, zet
$a_{0} = \sqrt{2}$, $b_{0} = 0$, $p_{0} = 2 + \sqrt{2}$ en itereer :
\\[ a_{k+1} = \frac{ \sqrt{a_{k}} + 1/ \sqrt{a_{k}} }{2} \]
\\[ b_{k+1} = \frac{ \sqrt{a_{k}} (1 + b_{k}) }{ a_{k} + b_{k} } \]
\\[ p_{k+1} = \frac{ p_{k} b_{k+1} (1 + a_{k+1}) }{ 1 + b_{k+1} } \]
\noindent
Dan convergeert $p_{k}$ kwadratisch naar $\pi$. De iteraties geven
vervolgens respectivelijk $3, 8, 19, 40, 83, 170, 345, 694, 1392,
2788 \ldots$ correcte decimalen.\\
\noindent
Voor het Borwein algorithme dat viermaal het aantal decimalen oplevert
per iteratie, start met $a_{0} = 6 - 4 \sqrt{2}$, $y_{0} = \sqrt{2} - 1$
en itereer dan :
\\[ y_{k+1} = \frac{ 1 - (1 - y_{k}^{4})^{1/4} } { 1 + (1 - y_{k}^{4})^{1/4} } \]
\\[ a_{k+1} = a_{k} (1 + y_{k+1})^{4} - 2^{2k+3} y_{k+1} (1 + y_{k+1} + y_{k+1}^{2}) \]
\noindent
Hiermee convergeert $a_{k}$ naar $\frac{1}{\pi}$.
\section{Andere Formules met $\pi$}
De twee natuurconstanten $e$ en $\pi$ zijn trouwens direct aan elkaar
gerelateerd. Een interessante relatie is :
\\[ e^{i \pi} + 1 = 0 \]
\noindent
Deze elegante formule met $0, 1, i$ en $\pi$ hebben we aan Euler te danken.
\noindent
Er zijn nog twee andere verassende formules die $e$ en $\pi$ met elkaar verbinden. Allereerst
\\[ e^{\pi} = 32 \prod_{j=0}^{\infty} \frac{a_{j+1}}{a_{j}} \]
\noindent
waarbij $a_{0}=\sqrt{2}$, $b_{0}=1$ en
\\[ a_{j}=\frac{a_{j-1} + b_{j-1}}{2} \]
\\[ b_{j}=\sqrt{a_{j-1} b_{j-1}} \]
\noindent
$a_{j}$ en $b_{j}$ worden respectievelijk arithmetisch en geometrisch
gemiddelde (AGM, Arithmatic Geometric Mean) genoemd. \\
\noindent
De volgende relatie is ook heel fraai :
\\[
\frac{e}{\pi} = \frac{ \sum_{j=0}^{\infty} \frac{1}{j!} }{ \sqrt{6\sum_{j=1}^{\infty}\frac{1}{j^{2}}} }
\]
\noindent
en volgt direct uit de Taylor expansie van $e$ in combinatie met
$\sum_{j=1}^{\infty} \frac{1}{j^{2}} = \frac{\pi^{2}}{6}$.
\section{Elementaire Arithmetische Operaties}
Zoals gezegd volgt uit de bovenstaande algorithmes dat je een aantal
elementaire arithmetische operaties zelf moet implementeren voor
{\bf arbitraire} precisie. De implementatie voor optellen en aftrekken
is slechts het zelf bij houden van een carry. Ik zal hier efficiente
algorithmes geven voor vermenigvuldigen, delen en tenslotte worteltrekken.
Een schat aan informatie over elementaire operaties is te vinden
in de eenmans encyclopedie 'The Art of Computer Programming' van
Donald Knuth~\cite{KNUTH} : een meesterwerk dat eigenlijk in geen
enkele boekenkast meer mag ontbreken.
\subsection{Vermenigvuldiging}
Jammer genoeg is de gewone 'lagere school methode' erg inefficient bij
een groot aantal decimalen, nl. orde $n^2$ voor het vermenigvuldigen van
twee getallen van $n$ decimalen.
De snelste~\cite{SCHOENHAGE,BAILEYFFT} vermenigvuldigingsmethoden maken
gebruik van F(ast) F(ourier) T(ransforms)~\cite{SWARZTRAUBER,SEDGEWICK}
maar ik zal simpelere alternatieven noemen. Ze zijn niet zo snel als FFTs,
maar wel redelijk eenvoudig te implementeren.
Als je maar een $n$-decimale nauwkeurigheid hebt en je moet de getallen
$A$ en $B$ vermenigvuldigen die beide $n$ decimalen groot kunnen zijn,
past het resultaat $C = A B$ niet altijd in $n$ decimalen, want dat kan
$2n$ decimalen groot zijn. In plaats van vermenigvuldigen van
afzonderlijke decimalen, kun je $A$ en $B$ opsplitsen in een $h$oog en
een $l$aag deel, zodat :
\begin{eqnarray}
C = & A B \nonumber \\
= & (10^{n/2} A_{h} + A_{l}) (10^{n/2} B_{h} + B_{l}) \nonumber \\
= & A_{h} B_{h} 10^{n} + ( A_{h} B_{h} + A_{l} B_{l} - (A_{h}-A{l})(B_{h}-B_{l}) ) 10^{n/2} + A_{l} B_{l} \nonumber
\end{eqnarray}
\noindent
Nu hoef je slechts de vermenigvuldigingen $A_{h} B_{h}$, $A_{l} B_{l}$ en
$(A_{h}-A{l})(B_{h}-B_{l})$ uit te voeren waarvan de uitkomsten allen
in $n$ decimalen passen. Vervolgens moet je schuiven en optellen.
Het bovenstaande algorithme kan recursief uitgevoerd worden door $A$
en $B$ herhaald in tweekn te splitsen totdat de stukken klein genoeg
zijn om met normale precisie de 3 vermenigvuldigingen uit te voeren.
Dit leidt tot een algorithme dat van de orde $n^{1.52}$ is.
Een andere mogelijkheid tot vermenigvuldigen wordt duidelijk wanneer
je beseft dat $A B = \frac{1}{2}(A+B)^{2} - A^{2} - B^{2})$.
Hiermee is ook een redelijke multi-precisie vermenigvuldiging te
implementeren omdat kwadrateren met reciproke operaties gedaan kan
worden wegens : $x^{2} = 1/(\frac{1}{x} - \frac{1}{x+1}) - x$
\subsection{Delen : de Reciproke}
Delen is een combinatie van een vermenigvuldiging met een reciproke
operatie. In plaats van $C = \frac{A}{B}$ bereken je dus eerst
$\frac{1}{B}$ en vermenigvuldigt vervolgens met $A$.
Vermenigvuldigen is hierboven al behandeld, dus hoe berekent men
$\frac{1}{B}$ ?
Aangezien $\frac{1}{B}$ het nulpunt is van de functie $f(x) = B -
\frac{1}{x}$, kun je het effectief met Newton iteratie berekenen :
het nulpunt van een functie $f(x)$ is de limiet $k \rightarrow \infty$
van $x_{k+1} = x_{k} - f(x_{k})/f'(x_{k})$.
Voor onze $f(x)$ moet je dus $x_{k+1} = x_{k} (2 - B x_{k})$ itereren.
Newton iteratie heeft het voordeel dat het aantal goede decimalen
verdubbelt bij iedere iteratie. Je hoeft dus niet gelijk vanaf het
begin van de berekening van $\frac{1}{B}$ met de volle nauwkeurigheid
te werken, maar verdubbelt deze bij iedere opvolgende iteratie.
\subsection{Worteltrekken}
Ook worteltrekken~\cite{ALT} laat zich het best door Newton iteratie temmen.
Het handigste is om niet $\sqrt{A}$ rechtstreeks met iteratie te berekenen,
maar $\frac{1}{\sqrt{A}}$ en vervolgens het geheel te completeren met
een laatste volle precisie vermenigvuldiging met $A$.
De Newton iteratie $x_{k+1} = \frac{x_{k} (3 - A x_{k}^{2})}{2}$ voor
$\frac{1}{\sqrt{A}}$ heeft namelijk het voordeel dat die geen deling meer
bevat (behalve delen door 2, maar dat is triviaal) en dat scheelt een stuk.
Natuurlijk laat dit alles zich generaliseren voor een $n$-de machts wortel
van $A$, zoek dan met Newton iteratie naar het nulpunt van
$f(x) = A - x^{n}$.
\section{Slot Opmerkingen}
In de praktijk kan er veel gewonnen worden door vanaf het begin decimaal
te rekenen met een eigen BCD-implementie omdat je uiteindelijk toch de
decimale representatie van $\pi$ wilt zien. Conversie van binair of
hexadecimaal naar decimaal is een zeer tijdrovende laatste operatie!
Wie niet zo geinteresseerd is in een zeer groot aantal decimalen of wie
gewoon een bloedsnelle computer heeft kan het onderstaande
C-programmaatje (voor 200 decimalen, maakt gebruik van de arctangens
reeks en getallenbases transformaties) gebruiken en/of aanpassen.
\begin{verbatim}
#include<stdio.h>
int b,c=2800,d,e=0,f[2801],g;
main()
{
while (b!=c) f[b++]=2000;
while (c!=0)
{
d=0 ; g=c*2-1; b=c;
while (b!=0)
{
d*=b; d+=f[b]*10000; f[b]=d%g; d/=g; g-=2; b--;
}
c-=14; printf("%.4d",e+d/10000); e=d%10000;
}
}
\end{verbatim}
\noindent
Voor de Borwein algorithmes is het natuurlijk handig als je van
tevoren beslag kunt leggen op de expansie van $\sqrt{2}$ met hetzelfde
aantal decimalen waarin je $\pi$ wilt berekenen, omdat $\sqrt{2}$ expliciet
verschijnt in de startwaardes van de iteraties~\cite{DUTKA}.\\
\noindent
Ter controle tenslotte, de eerste 100 decimalen van $\pi$ zijn :
\begin{verbatim}
3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510
58209 74944 59230 78164 06286 20899 86280 34825 34211 70679
\end{verbatim}
\noindent
De decimalen 990 t/m 1000 zijn {\tt 92164201989}
en de decimalen 4999990 t/m 5000000 zijn {\tt 20764619715}.
\begin{thebibliography}{99}
\bibitem{ALT}
H. Alt, ``Square Rooting Is as Difficult as Multiplication'', \\
Computing 21, pp. 221-232 (1979). \\
\bibitem{BAILEYFFT}
D.H. Bailey, ``A high performance fast Fourier transform algorithm for \\
the CRAY-2'', Journal of Supercomputing, v. 1, pp. 43-60 (1987). \\
\bibitem{BAILEYPI}
D.H. Bailey, ``The Computation of pi to 29,360,000 Decimal \\
Digits Using Borweins' Quartically Convergent Algorithm'', \\
Math. of Comp., v. 50, pp. 283-296 (1988). \\
\bibitem{BECKMANN}
P. Beckmann, ``A History of Pi'', Golem Press, Boulder, CO, 1971. \\
\bibitem{BORAGM}
J. M. Borwein, P.B. Borwein, ``The arithmatic-geometric mean and fast \\
computation of elementary functions'', SIAM Rev., v. 26, pp. 351-366 (1984). \\
\bibitem{BORPI}
J.M. Borwein, J.P. Borwein, ``More quadratically converging \\
algorithms for pi'',Math. Comp., v. 46, pp. 247-253 (1986). \\
\bibitem{BRENT}
R.P. Brent, ``Fast multiple-precision evaluation of elementary \\
functions'', J. Assoc. Comput. Mach., v. 23, pp. 242-251 (1976). \\
\bibitem{CIJSOUW}
P.L. Cijsouw, ``Commensurabiliteit en kettingbreuken'', Zenit, \\
Juni 1986, in Sterrenkunde op de huiscomputer, p. 213. \\
\bibitem{DUTKA}
J. Dutka, ``The Square Root of 2 to 1,000,000 Decimals'', \\
Math. of Comp., v. 25, pp. 927-930 (1971). \\
\bibitem{GAUSS}
C.F. Gauss, ``Werke'', Goettingen, 1863; 2nd ed., 1876, v. 2, p. 499-502. \\
\bibitem{KANADA}
Y. Kanada, Y. Tamura, ``Calculation of pi to 10,013,395 \\
Decimal Places Based on the Gauss-Legendre Algorithm and Gauss \\
Arctangent Relation'', Computer Centre, University of Tokyo, 1983. \\
\bibitem{KNUTH}
D. Knuth, ``The Art of Computer Programming'', Vol. 2: Semi- \\
numerical Algorithms, Addison-Wesley, Reading, Mass., 1981. \\
\bibitem{NEWMAN}
D. J. Newman, ``A Simplified Version of the Fast Algorithms of \\
Brent and Salamin'', Math. of Comp., v. 44, pp. 207-210 (1985). \\
\bibitem{RAMANUJAN}
S. Ramanujan, ``Modular equations and approximations to pi'', \\
Quart. J. Pure Appl. Math., v. 45, 1914, p.350-372; Collected papers \\
of Srinivasa Ramanujan, Cambridge 1927, p. 29-39. \\
\bibitem{RIETWEISNER}
G. Rietweisner, ``An ENIAC determination of pi and e to more \\
than 2000 decimal places", MTAC, v. 4, 1950, p. 11-15. \\
\bibitem{SAGAN}
Carl Sagan, ``Contact'', Arrow Books, 1986, ISBN 0-09-046950-2. \\
\bibitem{SALAMIN}
E. Salamin, ``Computation of pi using arithmetic-geometric mean'', \\
Math. of Comp., v. 30, pp. 565-570 (1976). \\
\bibitem{SCHOENHAGE}
A. Schoenhage, V. Strassen, ``Schnelle Multiplikation grosser \\
Zahlen'', Computing 7, pp. 281-292 (1971). \\
\bibitem{SEDGEWICK}
R. Sedgewick, ``Algorithms'', ISBN 0-201-06673-4, Chapter 41 : \\
The Fast Fourier Transform, pp. 583-593. \\
\bibitem{SHANKS}
D. Shanks, J.W. Wrench, Jr., ``Calculation of pi to 100,000 \\
decimals'', Math. of Comp., v. 16, pp. 76-99 (1962). \\
\bibitem{SWARZTRAUBER}
P. Swarztrauber, ``FFT algorithms for vector computers'', Parallel \\
Comput., v. 1, pp. 45-64 (1984).
\end{thebibliography}
\end{document}
Dit topic is gesloten.
:strip_exif()/u/22886/smileyoulooklikeshit.gif?f=community)
:strip_icc():strip_exif()/u/33320/crop5db40be9a5b1c_cropped.jpeg?f=community)
:strip_exif()/u/95433/shifty2.gif?f=community)
:strip_icc():strip_exif()/u/25814/camus_small.jpg?f=community)