Pi with Machin's formula (Erlang)

From LiteratePrograms

Jump to: navigation, search
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

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 regular integer.

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) ->
    Start = Unity div X,
    arccot(X, Unity, 0, Start, 1, 1).
arccot(_X, _Unity, Sum, XPower, N, _Sign) when XPower div N =:= 0 ->
    Sum;
arccot(X, Unity, Sum, XPower, N, Sign) ->
    Term = XPower div N,
    arccot(X, Unity, Sum + Sign*Term, XPower div (X*X), N+2, -Sign).

Applying Machin's formula

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

<<pi>>=
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>>=
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) ->
  Product;
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: 

<<machin.erl>>=
-module(machin).
-export([pi/1]).
arccot
pi
pow

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).
{ok,machin}
2> machin:pi(100).
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
Download code
Retrieved from "http://literateprograms.org/Pi_with_Machin%27s_formula_%28Erlang%29"
Views