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())
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