Recursion: Salamin and Brent equation for finding $\pi$

Question

I am attempting to compute $\pi$ using the Salamin and Brent equation:

• Salamin and Brent (1976): Set $a_0=1$, $b_0=\frac1{\sqrt{2}}$, $s_0=\frac12$. Then iterate: $$\small a_k=\frac{a_{k-1}+b_{k-1}}{2},\; b_k=\sqrt{a_{k-1} b_{k-1}},\; c_k=a_k^2-b_k^2,\; s_k=s_{k-1}-2^kc_k,\; p_k=\frac{2a_k^2}{s_k}$$ Then $p_k$ gives an approximation to $\pi$.

The thing is, this equation is very different from the ones I attempted before as it doesn't have a typical summation.

To start this, I attempted to define a[0] = 1 and b[0]=1/sqrt[2] etc. However, this doesn't seem right as the letters turn green and I'm not able to assign that particular value to the aforementioned variables.

However, since I'm new to this kind of computing, I am just unsure about how to approach computing a sum like this. What are some codes I need to learn in order to be able to do this?

Answer 1

You can implement recursions almost in the same way you have specified them:

Clear[a, b, c, s, p];
a[0] = 1.; b[0] = 1/Sqrt[2]; c[0] = 1/2; s[0] = c[0];
a[k_] := a[k] = (a[k - 1] + b[k - 1])/2;
b[k_] := b[k] = Sqrt[a[k - 1] b[k - 1]];
c[k_] := c[k] = a[k]^2 - b[k]^2;
s[k_] := s[k] = s[k - 1] - 2^k c[k];
p[k_] := p[k] = 2 a[k]^2/s[k];

The first 10 terms:

p[#] & /@ Range[10]

{3.18767, 3.14168, 3.14159, 3.14159, 3.14159, 3.14159, 3.14159, 3.14159, 3.14159, 3.14159}

It seems to converge very quickly! If you want to retain exact expressions (instead of machine precision), use a[0]=1 (instead of a[0]=1.). With this change you can get many decimal places of accuracy, for instance, to get the 10th iterate to 100 places:

N[p[10], 100]