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

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?

• sqrt or Sqrt? All built-in functions and symbols are capitalized. Commented Feb 20, 2021 at 22:49
• Commented Feb 20, 2021 at 22:51
• Crossposted here. Commented Feb 21, 2021 at 3:03

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]

• The code c[k_] := c[k] = c[k-1]^2/(16 a[k]^2); is better since it avoids loss of precision (at least for machine floating point). Somehow, Mathematica avoids this pitfall when using N[c[k],16] for example. Commented Feb 22, 2021 at 0:11
• As a matter of style, p[#] & /@ Range[10] is more naturally written as Array[p, 10]. Commented Feb 24, 2021 at 9:30