Ramanujan and Salamin-Brent Methods for π

MATH 242 Modern Computational Mathematics

Ramanujan's Method

Here is Ramanujan's sum (equation (1.12) in the text):
(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},
Do[
sum+=(((4k)!(1103+26390k))/((k!)^4*396^(4k))),
{k,0,n}];
sum*=2Sqrt[2]/9801;
1/sum (* return value *)
]

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

ramanujan[3]

N[ramanujan[3],50]

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

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

piValuesRamanujan[[1;;5]]

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

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

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.23 on page 36 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 π?

Notes about subscripted variables

Mathematica handles subscripts in a strange way. You can define subscripts like this:

(Subscript[a, 1])=7

...but the value isn’t associated with the variable that you expect.

?a

Instead, the value is stored as a “transformation rule” associated to the Subscript symbol:

DownValues[Subscript]

Instead of subscripts, use indexed variables, like this:

b[1]=7

b[2]=10

Now the values you stored are associated with the symbol that you expect:

?b

In your module, you only need to put the “head” of the indexed variable in the local variable list, like this: {a,b,s,p}. You can then assign indexed values such as a[0] and b[3] in your module.

Here is an implementation of the Salamin-Brent algorithm:

salaminBrent[m_]:=Module[
(* initialize variables *)
{a,b,s,p},
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,m}
];

(* return value *) 
p[m]
]

Compute the first two iterations:

salaminBrent/@Range[2]

Compute five iterations and convert the results to decimal values:

N[salaminBrent/@Range[5],20]

Compute the accuracy of the Salamin-Brent algorithm:

piValuesSalaminBrent=salaminBrent/@Range[8];

$MaxExtraPrecision=1000

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!