Ramanujan and Salamin-Brent Methods for π

Math 242 Modern Computational Mathematics

Ramanujan's Method

Here is Ramanujan's sum:
(1/π)=((2(Sqrt[2]))/9801)  (∑_(k=0)^∞)(((4k)!(1103+26390k))/((((k!))^4)(396^(4k))))

Complete the following module to add up the first n terms of Ramanujan’s sum and return an approximation of π:

ramanujan[n_]:=Module[
{sum=0},
(* main loop *)
Do[
sum+=(((4k)!(1103+26390k))/((k!)^4*396^(4k))),
{k,0,n}];
(* return value *)
(2*Sqrt[2]*sum/9801)^(-1)
]

Demonstrate that your module returns exact values that are close to π.

ramanujan[3]

Table[ramanujan[n],{n,1,4}]

Shorthand for Table is /@, like this:

ramanujan/@Range[4]

How many correct digits of π are produced by each term of Ramanujan’s sum?

nMax=20;
piValuesRamanujan=ramanujan/@Range[nMax];

$MaxExtraPrecision=1000;
accuracyRamanujan=N[-Log[10,Abs[Pi-piValuesRamanujan]],5]

ListPlot[accuracyRamanujan,AxesLabel->{"n","correct digits"}]

slope=(accuracyRamanujan[[20]]-accuracyRamanujan[[1]])/19

TableForm[Transpose[{Range[nMax],accuracyRamanujan}],TableHeadings->{None,{"terms","correct digits"}}]

Salamin-Brent Algorithm

Implement a module that computes approximations of π using the Salamin-Brent method (Algorithm 1.20 on page 31 in the text). Your module should accept one input parameter: an integer that specifies how many iterations of the algorithm to compute. Make a table showing how the number of correct digits of π increases with the number of terms that you add up. Answer the question: What is the relationship between the number of iterations and the number of correct digits of π?

A note about subscripts

Subscripts are good for presentation, not for variables!

(Subscript[a, 1])=3

?a

DownValues[Subscript]

Instead, use indexed variables:

a[1]=7

a[1]

a[2-1]

?a

Salamin-Brent Algorithm

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]
]

Compute the first two iterations:

salaminBrent/@Range[2]

Compute the accuracy:

piValuesSalaminBrent=salaminBrent/@Range[8];

accuracySalaminBrent=N[-Log[10,Abs[Pi-piValuesSalaminBrent]],5]

TableForm[Transpose[{Range[8],accuracySalaminBrent}],TableHeadings->{None,{"terms","correct digits"}}]

The Salamin-Brent algorithm approximately doubles the number of correct digits per iteration!