In-Class Demo

Math 242 Modern Computational Mathematics

Computing Digits of π

What does it mean to measure the time required to compute digits of π?

Consider the Salamin-Brent algorithm from last time. Here is an implementation:

salaminBrent[numIterations_]:=Module[
(* local variables *)
{a,b,s,p},

(* initial values *)
a[0]=1;
b[0]=1/Sqrt[2];
s[0]=1/2;

(* main loop *)
Do[
a[k]=(a[k-1]+b[k-1])/2;
b[k]=Sqrt[a[k-1]*b[k-1]];
s[k]=s[k-1]-2^k*(a[k]^2-b[k]^2);
p[k]=2a[k]^2/s[k],
{k,1,numIterations}
];

(* return value *)
p[numIterations]
]

Notice that the Salamin-Brent algorithm returns an algebraic expression involving fractions and exponents. It would take a lot more work to obtain actual digits from this expression.

salaminBrent[3]

Suppose that we want to compute 40,000 digits of π. First let’s find out how many iterations of the Salamin-Brent algorithm we will need for this. Compute a few iterations and look for a pattern:

piValuesSB=salaminBrent/@Range[12];
Block[{$MaxExtraPrecision=50000},
accuracySB=N[-Log[10,Abs[Pi-piValuesSB]],5]
]

Each iteration roughly doubles the number of correct digits. The 12th iteration gives us more than 5,500 correct digits, and we will need three more doublings after that. So we will need 15 iterations. Let’s see how long it takes to run the Salamin-Brent algorithm for 15 iterations:

ClearSystemCache[];
Timing[salaminBrent[15]][[1]]

Now let’s see how long it takes to run the Salamin-Brent algorithm for 15 iterations and also convert the output to a decimal number with 40,000 digits:

ClearSystemCache[];
Timing[N[salaminBrent[15],40000]][[1]]

We see that it takes much longer to compute the 40,000 digits than to compute 15 iterations of the Salamin-Brent algorithm!

"Dart Board" Algorithm for Approximating π

First we need to know how to generate random numbers.

RandomReal[]

Table[RandomReal[],20]

Now we can write a module that uses the dart board method to approximate π.

dartBoardPi[m_]:=Module[
(* local variables *)
{acc=0,x,y},

(* main loop *)
Do[
x=RandomReal[];
y=RandomReal[];
(*Print[{x,y}];*)

(* if (x,y) is within 1 unit of the origin, then increment acc *) 
If[x^2+y^2<1,acc+=1]
,m (* repeat the loop m times *)
];

(* return value *) 
4*acc/m
]

Try it out:

N[dartBoardPi[1000]]

Table[N[dartBoardPi[1000]],10]

Question 1:  If you use the dart board algorithm with m=1000, what percentage of the time does your approximation of π begin with the digits 3.1?

Here is one solution similar to what most people did in class:

proportion31[numTrials_]:=Module[
{acc=0, approx},

(* do the algorithm numTrials times,
count how many approximations are between 3.1 and 3.2 *)
Do[
approx=dartBoardPi[1000];
If[3.1<=approx<3.2,acc+=1],
numTrials
];

(* return the proportion *)
acc/numTrials
]

proportion31[100]

Now do 10,000 trials and see how many are between 3.1 and 3.2:

N[proportion31[10000]]

It appears that the approximation of π begins with the digits 3.1 about 66% of the time.

Bonus: Here is a fancy solution using the Select[] function

First use the dart board method to generate some approximations of π:

numTrials=10000;
approxPiVals=Table[dartBoardPi[1000],numTrials];

Now use Select[] to get only those approximations that are between 3.1 and 3.2. (For details, see the Mathematica documentation for Select.)

vals31=Select[approxPiVals,3.1<=#<3.2&];

Now compute the proportion:

N[Length[vals31]/numTrials]

We again found a percentage close to 66%.

Question 2: How many points (darts) are required to obtain an estimate of π that is within 0.1 of π at least 75% of the time?

One solution is to use trial and error. A module that performs a number if trials and returns the proportion of the approximations that are within 0.1 or[ will help:

proportionTenth[numTrials_, numDarts_]:=Module[
{acc=0, approx},

(* do the algorithm numTrials times with numDarts darts,
count how many approximations are within 0.1 of π *)
Do[
approx=dartBoardPi[numDarts];
If[Abs[Pi-approx]<0.1,acc+=1],
numTrials
];

(* return the proportion *)
acc/numTrials
]

Now try it for various values of darts with (say) 1000 trials each:

N[proportionTenth[1000,100]]

N[proportionTenth[1000,200]]

N[proportionTenth[1000,300]]

N[proportionTenth[1000,400]]

It appears that the number of darts is between 300 and 400. We could narrow this down a bit more, this time using a Table. Using a large number of trials will reduce the random variation in the results.

vals=Table[{n,N[proportionTenth[100000,n]]},{n,300,400,10}]

TableForm[vals,TableHeadings->{{},{"n","proportion"}}]

It looks like the number of darts required is somewhere around 360 to 380.