hoe werkt de BBP-formule voor Pi te berekenen?

manu, 13 jaar
13 juni 2010

ik ben al een tijdje op zoek naar een formule om pi te berekenen met <1000 getallen na de komma.
ik heb een script gemaakt in php met de BBP-formule.
$repeat=0;
$k=0;
$pi=0;
while($repeat!=10){
$iteration=(1/pow(16,$k))*((4/(8*$k+1))-(2/(8*$k+4))-(1/(8*$k+5))-(1/(8*$k+6)));
$pi=$pi+$iteration;
$repeat++;
$k++;
}
echo $pi;
?>
dit geeft als uitkomst : 3,1415926535898
als ik meer iteraties toevoeg of niet dan blijft het hetzelfde.
hoe doe je dat dan met 100 getallen na de komma ?
daarom dacht ik aan het BBP-algoritme, ze zeggen dat dat het n'de getal na de komma in het 16-tallige talstelsel kan geven, maar dat lukt mij niet.(ik zou willen dat het script de gatallen opslatt in een .txt bestandje)
wat doe ik fout ?

Antwoord

Hallo Manu,

Sorry dat het een paar dagen tijd heeft gevergd, maar je krijgt dan ook een uitgebreid antwoord van me.

1) Laat me eerst uitleggen waarom je met de code die je geeft niet meer dan 13 decimalen van pi kan berekenen.

Om dat te begrijpen moet je begrijpen hoe een computer met getallen omgaat. De getallen waarmee je werkt worden tijdens de berekeningen (optellen, aftrekken, vermenigvuldigen, delen) opgeslagen in het RAM (= random access memory) geheugen van de computer; in het geheugen worden vakjes gemaakt, typisch van 32 of 64 bits (dat betekent dat je in één "vakje" 32 of 64 nullen of enen kan opslaan). Per getal kan je dus maar een bepaalde hoeveelheid informatie, dus een bepaald aantal decimalen opslaan.

Meestal gaat men getallen opslaan onder de vorm van zogenaamde "floating point" getallen (zoek zelf eens "zwevendekommagetal" op in Wikipedia). Dit betekent dat een getal n wordt geschreven in "wetenschappelijke notatie" als

n = m x 10e

m is de mantisse en e is de exponent. In een typische floating point van het type "double" worden er 64 bits gebruikt (64 nullen of enen); de mantisse wordt opgeslagen met 53 bits en de exponent met 11 bits. En nu komt het: als er maximum 53 bits (in het binaire stelsel) zijn voor de mantisse, dan kan de mantisse in het decimale stelsel maar ongeveer 16 beduidende cijfers hebben.

Jij vindt pi als 3,1415926535898 met 14 beduidende cijfers. De laatste 2 cijfers laat de computer weg omdat er bij het rekenen ook kleine afrondingsfouten optreden.

Als je dus meer termen uit de BBP-reeks berekent, dan gebeurt er zoiets als hetvolgende

pi = 3,1415926535898 - 0,00000000000001

Op dat moment kan de computer de som niet meer verder uitwerken omdat er niet genoeg plaats is voor meer decimalen in het floating point getal.

Dat is waar het in jouw code misloopt. Is dat duidelijk?

2) Toch kunnen wetenschappers pi met veel meer decimalen berekenen. Bijvoorbeeld pi met 100 decimalen is

pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170680

Hoe kan je dit nu doen? Eén manier is om heel speciale floating point getallen te gebruiken met veel meer bits voor de mantisse. Dit kan je echter niet doen met de eenvoudigste programmeertalen en zulke grote floating point getallen vergen heel complexe code.

Een andere manier, die pas recent werd gevonden, is met de BBP-formule waarmee je een bepaalde decimaal van pi kan berekenen.

3) Laten we proberen te begrijpen hoe de BBP-formule werkt. Eerst en vooral geeft de BBP-formule de n-de digit van pi in het hexadecimale stelsel, ik heb het even voor jou berekend en pi in het hexadecimale stelsel is (eerste 17 digits):

pi = 3,243f6a8885a308d3

Stel nu dat we de 2de digit van pi willen berekenen (dat is dus "4"). Hiervoor vermenigvuldigen we de BBP-formule met 16. De redenering hier achter is dat we nu

16 x pi = 32,43f6a8885a308d3    (HEXADECIMAAL)

hebben; door te vermenigvuldigen met 16 hebben we de digits naar links verschoven zodat de digit die we zoeken nu net rechts van het kommateken staat.

Vervolgens berekenen we de eerste 2 termen van de BBP-formule:

16 x pi = 16 x (4 - 1/2 - 1/5 - 1/6) + (4/9 - 2/12 - 1/13 - 1/14) = 50,2628      (DECIMAAL)

We verwijderen nu het gehele deel van het resultaat:

50,2628 -> 0,2628

en we vermenigvuldigen nog eens met 16 om de digit die we zoeken voor het kommateken te krijgen:

0,2628 x 16 = 4,2208

Wat nu voor het kommateken ("4" ) staat is de 2de digit van pi die we zochten!

Om de derde digit te berekenen moet je pi vermenigvuldigen met 16^2 (om de gezochte digit net achter de komma te brengen), de eerste 3 termen van de BBP-formule uitrekenen, en tenslotte nog eens met 16 vermenigvuldigen om de gezochte digit voor de komma te brengen.

Om de vierde digit te berekenen moet je pi vermenigvuldigen met 16^3 (om de gezochte digit net achter de komma te brengen), de eerste 4 termen van de BBP-formule uitrekenen, en tenslotte nog eens met 16 vermenigvuldigen om de gezochte digit voor de komma te brengen.

Begrijp je het systeem? Probeer zelf eens (met de hand, nog niet met de computer!) de derde en de vierde digit te bepalen aan de hand van de BBP-formule om te tonen dat je de redenering begrijpt.

4) Dit lost echter nog niet het floating point probleem op, want je zal natuurlijk hetzelfde probleem hebben bij het optellen van de verschillende termen van de BBP-formule als je pakweg de 100ste hexadecimale digit van pi wil berekenen.

De truc om het floating point probleem te vermijden is in de som

16 x (4 - 1/2 - 1/5 - 1/6) + (4/9 - 2/12 - 1/13 - 1/14) = 50,2628 -> 0,2628

bij elke term eerst het gehele deel weg te gooien, dus wat je moet doen voor de tweede digit is:

      16 x (4   - 1/2 - 1/5 - 1/6          ) + (4/9           - 2/12         - 1/13         - 1/14         )
->          (64 - 8    - 3.2 - 2.6666...) + (0.4444... - 0.1666... - 0.0769... - 0.0714... )
->          (0   - 0    - 0.2 - 0.6666...) + (0.4444... - 0.1666... - 0.0769... - 0.0714... )

Je telt nu de termen van links naar rechts op, maar nadat je 2 termen hebt opgeteld, gooi je het gehele deel opnieuw weg of tel je er een geheel getal bij op om het resultaat terug positief te maken, dus

0 - 0 - 0.2 -> -0.2 -> 0.8

0.8 - 0.6666 = 0.1334 -> 0.1334

0.1334 + 0.4444 = 0.5778 -> 0.5778

enzoverder, uiteindelijk vind je 0.2629. Als je dit x 16 doet vind je 4.2064, en het gehele deel is de digit die je zocht.

5) Tenslotte moeten we een manier hebben om het gehele deel weg te gooien, zonder eerst het volledige getal te berekenen! Anders blijft hetzelfde probleem van de floating-point getallen zich stellen bij het berekenen van elke term. Gelukkig is er een algoritme gekend om de rest bij deling te berekenen voor gehele getallen. Dus als

a = b^c

met b en c allebei gehele getallen, dan kan je de rest bij deling door een ander geheel getal m berekenen zonder eerst b^c te berekenen. Ik noteer de rest bij deling als (a mod m)

Voorbeeld: (2^4 mod 3) = (16 mod 3) = 1 want 16 = 5x3 + 1

De 1-ste term van de BBP-formule (na vermenigvuldiging met 16^(n-1) als je de n-de digit wil bepalen, ziet eruit als volgt:

16^(n-1) / (8k + 1)

We kunnen het gehele deel verwijderen door in de teller de rest bij deling door de noemer te berekenen, dus wat we eigenlijk moeten hebben is

( 16^(n-1) mod (8k + 1) ) / (8k + 1)

( 16^(n-1) mod (8k + 1) ) is de rest bij deling van het geheel getal 16^(n-1) door het geheel getal 8k + 1. Zoals ik zei hebben we daar een truukje voor. Je kan hetzelfde resultaat bekomen door hetvolgende te berekenen:

16 x 16 x 16 x ... x 16     ( n - 1 factoren )

en telkens je x 16 doet neem je de rest bij deling door 8k + 1

Bijvoorbeeld (n = 4 en k = 1):

    16^3 mod 9 = ( (16 mod 9) x 16 x 16 ) = (7 x 16 x 16) mod 9 = ( (112 mod 9) x 16 ) mod 9
= (4 x 16) mod 9 = 64 mod 9 = 1

controle: 16^3 = 4096 = 455x9 + 1

Iets gelijkaardig kan je doen voor alle andere termen in de BBP-formule.

6) Je ziet dat het niet zo voor de hand liggend is. Ik weet dat er heel veel informatie in deze tekst staat. Print alles uit en neem je tijd om de hele tekst te overlopen; probeer elke stap te begrijpen en reageer als je een van de stappen echt niet begrijpt. Eens je alle stappen begrijpt, dan kan je proberen om een computercode te schrijven om de 100ste digit van pi te berekenen (in het hexadecimale stelsel). Je kan reageren op dit antwoord als je het resultaat hebt en dan zal ik zeggen of je juist bent. :-) Veel geluk!

Groetjes,

Philippe

Onderzoeker Fysica, Vrije Universiteit Brussel & Iowa State University

Reacties op dit antwoord

Er zijn nog geen reacties op deze vraag.

Enkel de vraagsteller en de wetenschapper kunnen reageren op een antwoord.

Beantwoord door

Prof. dr. Philippe Tassin

toegepast natuurkunde; optica; fotonica; fysica

Vrije Universiteit Brussel
Pleinlaan 2 1050 Elsene
http://www.vub.ac.be/

Zoek andere vragen

© 2008-2022
Ik heb een vraag wordt gecoördineerd door EOS vzw