My inquiry into a formula for calculating PI was piqued by the
Scientific American article in 1988 by Borwein, about, "Ramanujan &
PI."  Borwein is one of two Canadian brothers (Math Dep't at Dalhousie
University, Halifax, Nova Scotia) who together supplied the
Quartically Convergent Algorithm that allowed David H. Bailey (of the
NASA Ames Research Centre) to calculate PI to 29,360,000 decimal places
in January of 1986.

It is one thing to come up with an algorithm that will converge
quickly on PI.  It is quite another to run that algorithm efficiently
on a computer - even the Cray supercomputer Bailey used.  There are a
few arithmetic problems to solve in one's program, before actually
running the thing.

So it was that I went to the Math Library at the University of
British Columbia in the summer of '93 on a search.  It was for the
Holy Grail of number theory.  PI.

Since it had long since been proved, by the "Gregory Series", that:

   ARCTAN(x) = x - (x^3)/3 + (x^5)/5 - (x^7)/7 ......

Then it is fairly easy to number-crunch either the John Machin (1706)
formula for PI which goes:

  PI/4 = 4(ARCTAN(1/5)) - ARCTAN(1/239)

Or the more rapidly converging Stormer variation:

  PI = 24(ARCTAN(1/8)) + 8(ARCTAN(1/57)) + 4(ARCTAN(1/239))

Using the Machin series on Jul 29/61, Shanks and Wrench calculated PI
to 100,000 places, using an IBM 7090, in 8 hours and 43 minutes.  They
said that to reach 1,000,000 digits using the Stormer algorithm and the
7090, it would take 10 months.

Well, computers have changed in the last 32 years!  So have the
mathematics involved....

At the Math library, I ran into one of the professors (about 15 years
younger than me!).  He said, "so you're a PI-hack, eh?  NOBODY uses the
arctan method anymore!"  He took me up to his office and gave me a copy
of David H. Bailey's Jan 10/86 paper describing how Bailey had used
Borweins' Quartically Convergent algorithms to calculate PI to
29,360,000 decimal places, using a NASA Ames Research Center Cray-2
supercomputer.  It took 40 hours.

Basically, the following is the 1984 algorithm by Borwein & Borwein
(used as a 'check' for the better Borwein algorithm):

Let a(0)=the square root of 2
Let b(0)=0
Let p(0)= 2 + the square root of 2

Then iterate the following formulae:

a(n+1) =    (square root (a(n)) + ( 1 / square root (a(n)))
             ---------------------------------------------
                                    2

b(n+1) =    square root (a(n)) * ( 1 + b(n))
             ------------------------------
                        a(n) + b(n)

p(n+1) =    p(n) * b(n+1) * ( 1 + a(n+1))
             -------------------------
                      1 + b(n+1)


p(n) converges quartically to PI, successive iterations of this
algorithm yield 3, 8, 19, 41, 83, 170, 345, 694, 1392, 2788....
correct digits.  This algorithm is NOT self-correcting, EACH of the
ten iterations (to produce 2788 'correct' digits) must be performed
with more than 2788 digits of precision.


In the other, and better, Borwein algorithm used by Bailey, you find
the inverse of PI (ie.  1/PI ).  Using Newton's formula for finding
reciprocals and reciprocal of square roots;

x(n+1) = 2*x(n) - A*(x(n)^2)   converges to the reciprocal of A

y(n+1) = 3*y(n) - A*(y(n)^3)
          -----------------
                  2            converges to the reciprocal of the
                                       square root of A



both of which are necessary for 'flipping' a 29,360,000 digit number
(!), you use the other Borwein (and best) convergence as follows:

      let   a(0) = 6 - 4*(2^(1/2))
            y(0) = 2^(1/2) - 1       (how do you do a 'square root'
                                      sign in ASCII?)

then   y(n+1) = (1 - (1-y^4)^(1/4)) / (1 + (1-y^4)^(1/4))

   a(n+1) = a(n)(1+y(n+1)^4 - ( 2^(2n+3) * y(n+1)(1+y(n+1)+y(n+1)^2))

a(n) converges to 1/PI quadratically, and after you 'flip' it using
Newton's reciprocal formula, the very first term delivers 12 decimal
digits for PI, the second 32, third 76, fourth 164, fifth 332, sixth
580, all the way to the tenth term delivering 11,152 decimal digits. 
It is amazing how fast it converges!

Bailey calls the array method of keeping track of multiplications:

z(0) = x(0)*y(0)
z(1) = x(1)*y(0) + x(0)*y(1)
z(2) = x(2)*y(0) + x(1)*y(1) + x(0)*y(2)

(....... all the way to .............)

z(n-2) = x(n/2 - 1) * y(n/2 - 1)
z(n-1) = 0

Then you release the 'carries', so that if z(m)>9,

then z(m-1) = z(m-1) + INT(z(m)/10)
     z(m)   = z(m)   - 10*INT(z(m)/10)


the "schoolboy" approach, as he has used multiprecision arithmetic to
reduce the number of CPU calculations by a factor of 5, while
eliminating the rounding errors I got with big numbers, and the
computer's inability to keep intetegers greater than 8 digits. 
Multiprecision arithmetic is *the* way to go - just as important as
finding the quickest algorithm.



BTW - I found Ramanujan's formula, not mentioned in Borweins'
Scientific American article, which has been surpassed only by the
quadratic convergency... his was:

1/PI = 1/4 * ( 1123/882 - (22583/882^3 * 1/2 * 3/4^2) + (.....)
It doesn't look that exciting, and the reason they don't use it (even
though it is slightly faster than Stormer's) is that the terms are
impossible to double check (such checking using a secondary formula
required to claim a legitimate result).  Ramanujan had 14 such
convergences.  Where he got them from, the best and brightest
mathematicians of today do not know.

Now, onto the Chudnovsky Bros. of Columbia University, the current
record holders at 2.26 billion digits.  Who is going to pry their
algorithm away from them?  It promises to allow many users to
participate in a PI 'chain letter', allowing an indefinite
calculation.  When one calculator gets bored, he or she simply passes
the whole thing to someone else.  Trillions and ga-zillions are
possible.  Boo-ha-ha-ha-ha-ha....


Stuart A. Lyster
