Pi with Machin's formula (Python)

Other implementations: Erlang | Haskell | Java | Lisp | Python

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.

Contents

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

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.

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


Applying Machin's formula

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

<<pi.py>>=
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)
31415926535897932384626433832795028841971693993751058209749445923078164062862089
986280348253421170679L


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

hijacker
hijacker
hijacker
hijacker