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?


1 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]
  • 1
    $\begingroup$ 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. $\endgroup$
    – Somos
    Commented Feb 22, 2021 at 0:11
  • $\begingroup$ As a matter of style, p[#] & /@ Range[10] is more naturally written as Array[p, 10]. $\endgroup$
    – Roman
    Commented Feb 24, 2021 at 9:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.