= Area: CALG. ECHO: Program ==================================================
  Msg#: 45                                           Date: 03 Nov 94  22:40:47
  From: Matthew Mastracci                            Read: Yes    Replied: No
    To: All                                          Mark:
  Subj: Super Pi!
==============================================================================
Hi All!

     Have you ever been stuck for an exact value of Pi?  Well, this program
can calculate Pi to roughly 2^16 digits in under a second!  (on a 286/10
without a coprocessor too!)  It converges to Pi exponentially, so it would
probably be easy to calculate Pi to 16 million digits (with adequate variable
and result storage, that is) in less than 5 mins.

{$N+,E+,R-}
program Pi;
{ Determines the exact value of Pi using a formula which converges it }
uses Crt;

var
  a, b, c : Extended;
  x, y : Extended;

begin
  a := 1;
  b := Sqrt(0.5);
  c := 1 / 4;
  x := 1;
  repeat
    y := a;
    a := (a + b) / 2;
    b := Sqrt(b * y);
    c := c - x * Sqr(a - y);
    x := 2 * x; { x is equal to the number of exact digits }
    GotoXY(1, 1);
    Write(Sqr(a + b) / (4 * c):20:20); { Converges to Pi }
    ClrEOL;
    WriteLn;
  until x >= 128; { Roughly 128 digit precision: overkill }
end.

Background info of the function: (based on a short article by David Wells)

     The most recent method for calculating Pi, which was used by Tamura and
Kanada for their calculation to 16 million places, is based on Gauss's study of
the arithmetic-geometric mean of two numbers.

     Instead of using an infinite sum or product, the calculation goes round
and round in a loop.  It has the amazing property that the number of correct
digits approximately doubles with each circuit of the loop, so that going round
a mere 19 times gives Pi correct to over 1 million decimal places!

     The steps must be followed in sequence, up to (A + B) ^ 2 / 4C, which is
the first approximation to PI.  Then return to the first step (Y = A) and go
around again.  Here is the simple loop for calculating Pi:

     Y = A
     A = (A + B) / 2
     B = SQRT(BY)
     C = C - X(A - Y)^2
     X = 2X
     Print (PI) = (A + B) ^ 2 / 4C

The initial values are, A = X = 1, B = 1 / SQRT(2), and C = 1 / 4.

Regards,
  /\/\att /\/\astracci
  rayban@matrix.cpubbs.cuug.ab.ca

-!- GoldED
 ! Origin: Matrix (1:134/135.905)