Pi with Machin's formula (Haskell)

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.

Contents


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

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

<<arccot>>=
arccot x unity =
    arccot' x unity 0 start 1 1
      where start = unity `div` x
            arccot' x unity sum xpower n sign | xpower `div` n == 0 = sum
                                              | otherwise           =
                arccot' x unity (sum + sign*term) (xpower `div` (x*x)) (n+2) (-sign)
                  where term = xpower `div` n

[edit] Applying Machin's formula

Finally, the main function, which uses Machin's formula to compute π using the necessary level of precision (the name "pi" conflicts with the pre-defined value "pi" in Prelude):

<<pi>>=
machin_pi digits = pi' `div` (10 ^ 10)
                     where unity = 10 ^ (digits+10)
                           pi' = 4 * (4 * arccot 5 unity - arccot 239 unity)

Now we put it all together in a module:

<<machin.hs>>=
arccot
pi

[edit] Running

$ ghci
GHCi, version 6.8.2: http://www.haskell.org/ghc/  :? for help
Loading package base ... linking ... done.
Prelude> :l machin.hs
[1 of 1] Compiling Main             ( machin.hs, interpreted )
Ok, modules loaded: Main.
*Main> machin_pi 100
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
Download code
hijacker
hijacker
hijacker
hijacker