Pi with Machin's formula (Python)

From LiteratePrograms
Jump to: navigation, search
Other implementations: Erlang | Haskell | Java | Lisp | Python

[edit] Machin's formula

A simple way to compute the mathematical constant π ≈ 3.14159 with any desired precision is due to John Machin. In 1706, he found the formula

\frac{\pi}{4} = 4 \, \arccot \,5 - \arccot \,239

which he used along with the Taylor series expansion of the arc cotangent function,

\arccot x = \frac{1}{x} - \frac{1}{3 x^3} + \frac{1}{5 x^5} - \frac{1}{7 x^7} + \dots

to calculate 100 decimals by hand. The formula is well suited for computer implementation, both to compute π with little coding effort (adding up the series term by term) and using more advanced strategies (such as binary splitting) for better speed.


In order to obtain n digits, we will use fixed-point arithmetic to compute π × 10n as a Python long.

[edit] High-precision arccot computation

To calculate arccot of an argument x, we start by dividing the number 1 (represented by 10n, which we provide as the argument unity) by x to obtain the first term. We then repeatedly divide by x2 and a counter value that runs over 3, 5, 7, ..., to obtain each next term. The summation is stopped at the first zero term, which in this fixed-point representation corresponds to a real value less than 10-n.

def arccot(x, unity):
    sum = xpower = unity // x
    n = 3
    sign = -1
    while 1:
        xpower = xpower // (x*x)
        term = xpower // n
        if not term:
        sum += sign * term
        sign = -sign
        n += 2
    return sum

[edit] Applying Machin's formula

Finally, the main function, which uses Machin's formula to compute π using the necessary level of precision:

def pi(digits):
    unity = 10**(digits + 10)
    pi = 4 * (4*arccot(5, unity) - arccot(239, unity))
    return pi // 10**10

To avoid rounding errors in the result, we use 10 guard digits internally during the calculation. We may now reproduce Machin's result:

>>> pi(100)

The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer. (More sophisticated techniques are necessary to calculate millions or more digits in reasonable time, although in principle this program will also work.)

Download code