Pi with Machin's formula (Python)
From LiteratePrograms
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
which he used along with the Taylor series expansion of the arc cotangent function,
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.)
Download code |