Pi with Machin's formula (Erlang)

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 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(X, Unity) ->
    Start = Unity div X,
    arccot(X, Unity, 0, Start, 1, 1).

arccot(_X, _Unity, Sum, XPower, N, _Sign) when XPower div N =:= 0 ->
arccot(X, Unity, Sum, XPower, N, Sign) ->
    Term = XPower div N,
    arccot(X, Unity, Sum + Sign*Term, XPower div (X*X), N+2, -Sign).

[edit] Applying Machin's formula

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

pi(Digits) ->
    Unity = pow(10, Digits+10),
    Pi = 4 * (4 * arccot(5, Unity) - arccot(239, Unity)),
    Pi div pow(10,10).

Erlang's math:pow(X,N) returns a float so it will not work for this calculation. Instead we implement a simple pow function which uses exponentiation by squaring:

pow(Base, Exponent) when Exponent < 0 ->
  pow(1 / Base, -Exponent);

pow(Base, Exponent) when is_integer(Exponent) ->
  pow(Exponent, 1, Base).

pow(0, Product, _Modifier) ->

pow(Exponent, Product, Modifier) when Exponent rem 2 =:= 1 ->
  pow(Exponent div 2, Product * Modifier, Modifier * Modifier);

pow(Exponent, Product, Modifier) ->
  pow(Exponent div 2, Product, Modifier * Modifier).

Now we put it all together in a module:


[edit] Running

$ erl
Erlang (BEAM) emulator version 5.5.2 [source] [async-threads:0] [kernel-poll:false]

Eshell V5.5.2  (abort with ^G)
1> c(machin).
2> machin:pi(100).
Download code