Ramanujan and Salamin-Brent Methods for π

MATH 242 Modern Computational Mathematics

Ramanujan's Method

Here is Ramanujan's sum (equation (1.12) on page 33 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 m terms of Ramanujan’s sum and return an approximation of π:

ramanujan[m_]:=Module[{sum=0},
Do[
sum+=(((4k)!(1103+26390k))/((k!)^4*396^(4k))),
{k,0,m}];
sum*=2Sqrt[2]/9801;
1/sum
]

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

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

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

$MaxExtraPrecision=500;

accuracyRamanujan=N[-Log[10,Abs[Pi-piValuesRamanujan]],5]

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

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

It looks like each additional term produces 8 more correct digits of π!

Indexed Variables in Mathematica

Mathematica lets you type variables with subscripts, but these are used for display purposes and not for creating a list of variables in a computation. For example, consider the following:

Clear[a];
(Subscript[a, 1])=3
(Subscript[a, 2])=7

(Subscript[a, 1])

We might that the previous code has defined two variables associated with the symbol a, but this is not the case. We can see this if we query the symbol a, like this:

?nMax

?a

If there were values assigned to the symbol a, we would see them in the previous output, but we don’t see any values there.

Instead, our subscripted variables are actually stored as patterns associated with the Subscript object, as we can see by running the following code.

DownValues[Subscript]

If you want to create indexed variables for use in a computation, you should put the index in square brackets following the variable name. (This is the same way that you would define individual values of a function.) Here is an example:

Clear[b]
b[1]=3
b[2]=7

Now these values are actually associated with the symbol b, which we can see as follows:

?b

We can also access these values by using an index variable within the square brackets, provided that this index variable evaluates to one of the indexes that we defined:

Do[
Print[  b[i]  ],
{i,1,2}]

If you want to define indexed variables in a module, only the “head” of the variables (that is, the name that occurs before the square brackets) needs to go in the local variable list. Here is an example:

demoModule[]:=Module[
{a},

(* define some indexed variables *)
a[1]=2;
a[2]=5;
a[4]=8;

(* return the values as a list *)
{a[1],a[2],a[4]}
]

demoModule[]

Let's clear the symbols we’ve defined in this demonstration:

Clear[a, b,demoModule, Subscript]

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 π?

Here is an implementation using indexed variables:

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

(* initializations *)
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 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!