{
Copyright (C) 1998-99 Free Software Foundation, Inc.

Author: Frank Heckenbach <frank@pascal.gnu.de>

Demo program for the GMP unit

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 2.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA.
}

program Pi;

uses GPC, GMP;

var
  n : MedCard;
  f : mpf_t;
  s : CString;
  e : mp_exp_t;

procedure mpf_arctan (var c : mpf_t; const x : mpf_t);
(* $$\arctan x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{2n+1}$$ *)
var
  p, mx2, c0, s : mpf_t;
  Precision, n : MedCard;
begin
  mpf_init_set (p, x);
  mpf_init (mx2);
  mpf_mul (mx2, x, x);
  mpf_neg (mx2, mx2);
  mpf_init (c0);
  mpf_init (s);
  mpf_set (c, x);
  Precision := mpf_get_prec (c);
  n := 1;
  repeat
    mpf_mul (p, p, mx2);
    mpf_div_ui (s, p, 2 * n + 1);
    mpf_set (c0, c);
    mpf_add (c, c, s);
    Inc (n)
  until mpf_eq (c0, c, Precision) <> 0;
  mpf_clear (s);
  mpf_clear (c0);
  mpf_clear (mx2);
  mpf_clear (p)
end;

procedure mpf_pi (var c : mpf_t);
(* $$4 \arctan \frac{1}{5} - \arctan \frac{1}{239} = \frac{\pi}{4}$$*)
var b : mpf_t;
begin
  mpf_set_ui (c, 1);
  mpf_div_ui (c, c, 5);
  mpf_arctan (c, c);
  mpf_mul_ui (c, c, 4);
  mpf_init_set_ui (b, 1);
  mpf_div_ui (b, b, 239);
  mpf_arctan (b, b);
  mpf_sub (c, c, b);
  mpf_mul_ui (c, c, 4);
  mpf_clear (b)
end;

begin
  case ParamCount of
    0 : begin
          Writeln ('Pi computation.');
          Write ('Enter the number of digits: ');
          Readln (n);
        end;
    1 : ReadStr (ParamStr (1), n);
    else
      Writeln (StdErr, 'Pi computation.');
      Writeln (StdErr, 'Usage: ', ParamStr (0), ' [number_of_digits]');
      Halt (1)
  end;
  mpf_set_default_prec (Round (n * Ln (10) / Ln (2)) + 16);
  mpf_init (f);
  mpf_pi (f);
  s := mpf_get_str (nil, e, 10, n, f);
  {$X+} Writeln ('Pi ~ ', PChars (s)^[1 .. e], '.', s + e); {$X-}
  Dispose (s);
  mpf_clear (f)
end.
