Mean π: Archimedes vs. Gauss - π computation through generalized means

Question

I've written this simplified code to compute Pi for educational/demonstration purposes.

These methods are based upon the generalized means: see a presentation on Pi and the AGM.

Archimedes' method gives linear convergence, which means you get two extra bits of precision per iteration.

Gauss's method gives quadratic convergence, which means your precision doubles every iteration.

One can wrap these methods in timers or print current results to see how they converge.

Archimedes' method could have been used to disprove a millennia of false claims about Pi. Alas, history.

import decimal


def pi_arc():
    """Archimedes c. ~230 B.C.E."""
    a, b = D(3).sqrt() / D(6), D(1) / D(3)
    pi = 0
    while True:
        an = (a + b) / 2
        b = (an * b).sqrt()
        a = an
        piold = pi
        pi = 2 / (a + b)
        if pi == piold:
            break
    return D(str(pi)[:-3])


def pi_agm():
    """Gauss AGM Method c. ~1800 A.D. """
    a, b, t = 1, D(0.5).sqrt(), 1 / D(2)
    p, pi, k = 2, 0, 0
    while True:
        an = (a + b) / 2
        b = (a * b).sqrt()
        t -= p * (a - an)**2
        a, p = an, 2**(k + 2)
        piold = pi
        pi = ((a + b)**2) / (2 * t)
        k += 1
        if pi == piold:
            break
    return D(str(pi)[:-3])


if __name__ == "__main__":
    prec = int(input('Precision for Pi: '))
    """Plus 3 for error"""
    decimal.getcontext().prec = prec + 3
    D = decimal.Decimal
    print(pi_arc())
    print(pi_agm())

Answer 1

Convergence testing

    if pi == piold:
        break

This is not usually done, because float equality has a lot of gotchas. In this case it's possible due to the numbers being Decimal, but if you need to move away from Decimal you're going to encounter issues.

Usually, convergence is measured as the absolute error decreasing below a chosen epsilon, a very small positive number. One advantage is that if you start checking for convergence to epsilon, your code will be compatible with arbitrary-precision math, where the two numbers will never equal each other exactly but you still need sane termination criteria.

Formatting/rounding

D(str(pi)[:-3])

This looks troublesome. You're converting a float to a string, and then selecting a certain number of fixed digits. Don't do this. Instead, just use the built-in round, which works with Decimals just fine.

Order of Operations

pi = ((a + b)**2) / (2 * t)

Exponentiation takes precedence over division, so you can drop the first pair of outer parens.

Import assignment

D = decimal.Decimal

Usually you shouldn't do this, and instead you should use normal import syntax:

from decimal import Decimal as D

Answer 2

A few more points:

• The precision of the decimal context is set at the top-level of the code, while the result of computation is rounded within the two methods. That is not good because in this case the methods actually have no knowledge about how the precision is adjusted based on user input. These two operations should be handled on the same level: either both on the top-level or both on the method level (the precision value needs to be passed down as a function argument). I personally prefer the latter option.

• Despite the precision of the decimals are increased from user input, for convergence testing that might not be necessary. For example, if a precision value of three is requested, checking abs(pi_old - pi) < Decimal("0.001") seems enough to me (since the \$\pi>1\$, only two digits after decimal points are needed).

• Convergence testing can be done more cleanly using math.isclose (thanks to @StevenJackson's comment)

• Keeping the full class name Decimal rather than shortening it to a single-letter name D improves code readability. Nowadays most IDEs can autocomplete long names therefore using short names does really save much time.

• Multiple classes / methods can be imported in the one statement.

  from decimal import Decimal, getcontext
  

• Since the code is for demonstration purpose, it might be better to keep it as close to the presented pseudocode / description as possible. For example, a more direct implementation of the core loop of the Gauss-Legendre algorithm on page 31 of the linked PDF could be like the following.

  from itertools import count
  from math import isclose
  
  a = Decimal(1)
  b = Decimal(0.5).sqrt()
  s = 1 / Decimal(4)
  pi = 0
  min_delta = Decimal("0.001")  # minimum difference between steps to continue loop
  
  for n in count():    # Use itertools.count to generate an infinite sequence starting from 0
      a_next = (a + b) / 2
      pi_old = pi
      pi = a_next ** 2 / s
      if isclose(pi, pi_old, rel_tol=0, abs_tol=min_delta):
          break
  
      b = (a * b).sqrt()
      s -= 2**n * (a - a_next)**2
      a = a_next