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

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())
• Is there a specific part of the code you'd like us to review? It's unclear if you'd like this to be reviewed or if you just want to post this for educational purposes. – Confettimaker Sep 11 at 0:46
• All? Open to all suggestions to optimize or make more presentable. – TheHoyt Sep 11 at 1:13
• Asking out of curiosity: what are a millenia of false claims about Pi? – lukeg Sep 11 at 16:55
• From the article: 'David Bailey has observed that there are at least eight recent papers in the “refereed” literature claiming that π = (14 −√2)/4 = 3.1464 · · · , and another three claiming that π = 17 − 8√3 = 3.1435 · · · .' – TheHoyt Sep 11 at 18:15
• @TheHoyt: That's really weird, since both are worse approximations than 22/7. Surely 22/7 should suffice for everyone's "wrong π" needs! – Nick Matteo Sep 11 at 22:38

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
• Excellent points, I will be implementing everything you have stated. The decimal equality is a luxury I'm taking here in Python to demonstrate a certain aspect of the AGM, but I agree with you. – TheHoyt Sep 11 at 3:33
• I would recommend docs.python.org/3/library/math.html#math.isclose specifically for convergence testing. It works for (mixtures of) all the types I've used in python, and relative tolerance is also a nice option. (especially if you're worried about things changing to floats at some point, isclose "just works" with any combination of floats and decimals as the values or as the tolerances) – Steven Jackson Sep 11 at 19:36

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

• I would always recommend docs.python.org/3/library/math.html#math.isclose for convergence testing. It handles comparisons even across types cleanly and there's not ever a compelling reason not to use it. – Steven Jackson Sep 11 at 19:38
• @StevenJackson Thanks. I've updated my answer. – GZ0 Sep 11 at 19:52