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. A simple way is to use RandomReal[].

RandomReal[]

Table[RandomReal[],10]

RandomReal[{-1,1},2]

We also need to know how to write conditional statements (also called if statements). Here is an example:

r=RandomReal[]
If[r<0.6,Print["Hello"]]

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[];
(* if (x,y) is within 1 unit of the origin, then increment acc *) 
If[x^2+y^2<1,acc+=1]
,m];

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

Try it out:

dartBoardPi[1000]//N

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 way to do it.

findPercentage[numTrials_]:=Module[
{acc=0},

(* Run the dartboard algorithm with 1000 darts many times. Keep track of how many times the result is between 3.1 and 3.2. *)
Do[
If[3.1<=dartBoardPi[1000]<3.2,acc+=1],
numTrials
];

(* return the proportion *) 
acc/numTrials
]

result=findPercentage[500]

N[result]

It looks like the dartboard algorithm with 1000 darts returns a value between 3.1 and 3.2 about 65% of the time.

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?

This is similar to Question 1, but now the number of points (darts) is not specified. It’s fine to experiment with different numbers of darts until you find a number of darts that gives an approximation satisfying |π-approx|<0.1 more than 75% of the time.