(*
  This program computes PI by the formula:
  pi / 4 = 4 * ArcTan(1/5) - ArcTan(1/239)

  Adapted from some program I found on the
  internet a long time ago that computed PI by:
  pi / 4 = ArcTan(1/2) + ArcTan(1/3)

  I have no idea who wrote the original program
  nor remember where I found it--there wasn't any
  credits or comments in the original C source.

  There are faster ways of computing PI, and
  this program probably could be optimized with
  all kinds of tricks, but it is sufficiently fast
  for computing PI using relativly portable and
  understandable routines that can be easily
  translated into other languages.

  Computing PI using FFTs are much much faster,
  but also much more complicated.  Using the
  ArcTan methods are fairly simple to implement,
  but not nearly as fast as the FFT methods.  Also,
  the pi / 4 = 4 * ArcTan(1/5) - ArcTan(1/239)
  formula is one of the faster of the various
  ArcTan formulas for computing PI.

  Author: M. Scott Reynolds
  Date: 09 September 2006
  Compiler: Free Pascal
*)

program pi;

uses
  SysUtils;

const
  (* Change this to compute more digits of PI *)
  PRECISION = 5000;

var
  i: LongInt;
  start, finish: Comp;

  (* These arrays can be ShortInt, but I found
     that LongInts resulted in an increase in speed. *)

  p: array [0..PRECISION] of LongInt;    (* p := arctan(1/5) * 4 *)
  q: array [0..PRECISION] of LongInt;    (* q := arctan(1/239) *)
  t: array [0..PRECISION] of LongInt;    (* workspace for ArcTanInv *)
  u: array [0..PRECISION] of LongInt;    (* workspace for ArcTanInv *)

(* p := s + t.

  Assumes low end of the array is 0 and high end
  is PRECISION.
*)
procedure Add(var p, s, t: array of LongInt);
var
  j: LongInt;
begin
  (*
    Note that this routine, under certain situations could
    end up in an array out of bounds error, due to j-1
    if j = 0.  However, that never happens when used in this
    program, thus eliminating the need for checking j before
    referencing a potential -1 in the array.  Used anywhere
    else, that may not be the case.
  *)

  for j := PRECISION downto 0 do
  begin
    p[j] := s[j] + t[j];
    if p[j] > 9 then
    begin
      Dec(p[j], 10);	(* Decrement by 10 *)
      Inc(p[j-1], 1);	(* Increment by 1 *)
    end
  end
end;

(* p := s - t.

  Assumes low end of the array is 0 and high end
  is PRECISION.
*)
procedure Subtract(var p, s, t: array of LongInt);
var
  j: LongInt;
begin
  (*
    Note that this routine, under certain situations
    could end up in an array out of bounds error, due to
    j-1 if j = 0.  However, that never happens when used
    in this program, thus eliminating the need for
    checking j before referencing a potential -1 in the
    array.  Used anywhere else, that may not be the case.
  *)

  for j := PRECISION downto 0 do
  begin
    p[j] := s[j] - t[j];
    if p[j] < 0 then
    begin
      Inc(p[j], 10);  (* Increment by 10 *)
      Dec(p[j-1], 1); (* Decrement by 1 *)
    end
  end
end;

(* p := s * multiplier.

   Assumes low end of the array is 0 and high
   end is PRECISION.
 *)
procedure Multiply(var p, s: array of LongInt; const multiplier: LongInt);
var
  i: LongInt;
  b: LongInt;
  carry: LongInt;
begin
  carry := 0;

  for i := PRECISION downto 0 do
  begin
    b := s[i] * multiplier + carry;
    p[i] := b mod 10;
    carry := b div 10;
  end;
end;

(* p := s / divisor.

   Assumes low end of the arrays is 0 and high
   end is PRECISION.
 *)
procedure Divide(var p, s: array of LongInt; const divisor: LongInt);
var
  i: LongInt;
  b: LongInt;
  remainder: LongInt;
begin
  remainder := 0;

  for i := 0 to PRECISION do
  begin
    b := 10 * remainder + s[i];
    p[i] := b div divisor;
    remainder := b mod divisor;
  end;
end;

(* p := s / divisor. Return true if p = 0.

   Assumes low end of the array is 0 and high
   end is PRECISION.
 *)
function DivideIsZero(var p, s: array of LongInt; const divisor: LongInt):Boolean;

var
  i: LongInt;
  b: LongInt;
  remainder: LongInt;
begin
  DivideIsZero := True;	(* Default condition *)
  remainder := 0;

  for i := 0 to PRECISION do
  begin
    b := 10 * remainder + s[i];
    p[i] := b div divisor;
    remainder := b mod divisor;
    if DivideIsZero then
      if p[i] <> 0 then
	DivideIsZero := False;
  end;
end;

(* p := ArcTanInv(s).

   Assumes low end of the array is 0 and high
   end is PRECISION.
 *)
procedure ArcTanInv(var p: array of LongInt; const s: LongInt);
var
  n: LongInt;
  ss: LongInt;
  isZero: Boolean;
begin
  ss := s * s;

  (* t, being a global array, was automatically initialized to zero. *)
  t[0] := 1;

  Divide(t, t, s);

  for n := 0 to PRECISION do
    p[n] := t[n];

  n := 1;
  isZero := false;      (* While t is not zero. *)
  while not isZero do
  begin
    Inc(n, 2);		 (* Increment by 2 *)
    Divide(t, t, ss);
    Divide(u, t, n);
    Subtract(p, p, u);

    Inc(n, 2);

    (* Divide and check for t = 0 *)
    isZero := DivideIsZero(t, t, ss);

    Divide(u, t, n);
    Add(p, p, u);
  end;
end;

(* Compute the value of PI and print. *)
begin
  start := TimeStampToMSecs(DateTimeToTimeStamp(Now));

  (* p :=  (ArcTanInv(5) * 4 - ArcTanInv(239)) * 4 *)
  ArcTanInv(p, 5);
  Multiply(p, p, 4);
  ArcTanInv(q, 239);
  Subtract(p, p, q);
  Multiply(p, p, 4);

  finish := TimeStampToMSecs(DateTimeToTimeStamp(Now));

  WriteLn('pi = ', p[0], '.');
  for i := 1 to PRECISION do
  begin
    Write(p[i]);

    (* Format output. *)
    if i mod 1000 = 0 then
    begin
      WriteLn;
      WriteLn;
    end
    else if i mod 50 = 0 then
      WriteLn
    else if i mod 10 = 0 then
      Write(' ');
  end;
  WriteLn;
  WriteLn(finish - start, ' ms to compute ',
      PRECISION, ' digits of PI.');
end.