NetNews Usenet Archive 1992 #23/NN_1992_23.iso/spool/sci/math/13299

Text File | 1992-10-16 | 11.6 KB | 237 lines

Newsgroups: sci.math
Path: sparky!uunet!cs.utexas.edu!sun-barr!ames!agate!stanford.edu!CSD-NewsHost.Stanford.EDU!Sunburn.Stanford.EDU!pratt
From: [e-mail] (Vaughan R. Pratt)
Subject: Re: Looking for fast methods of computing PI
Message-ID: <[e-mail]>
Sender: [e-mail]
Organization: Computer Science Department,  Stanford University.
References: <[e-mail]> <[e-mail]>
Date: Fri, 16 Oct 1992 11:12:22 GMT
Lines: 225

In article <[e-mail]> [e-mail] writes:
>I don't know what method you are talking about but I have been
>thinking this problem for a long time and I found a way(the hard way)
>to find PI.  I use the idea of polygon in a circle.  I first use
>a square to start with.  Then I turn the square into a octagon and find
>the length of the side in terms of square root of 2. And then I
>double the number of side again and again and again.....
>It works fine but the algebra is heavy.  I figure square root of 2
>is easier to find than PI.  A friend of mine used taylor's series
>to find Pi up to 5000 digit.

I'd like to take this opportunity to publicize some little-known
aspects of this method.  What is well-known is that this is Archimedes'
half-angle method, starting with a square rather than a hexagon.
Archimedes' method yields exactly one bit of precision per "major"
operation (= multiplication or square root), in contrast to Gene
Salamin's early-1970's application of the Lagrange-Gauss
arithmetic-geometric-mean method, which doubles the number of bits of
precision at each iteration.

Before Lagrange-Gauss-Salamin, all improvements over Archimedes were
mere constant factors, as follows.

    Inventor      Speed in bits/op     Method

        Archimedes, c.250 BC    1       half-angle method
    Vieta, 1593        1    product(r=sqrt((r+1)/2))
    Descartes, < 1650       1       constant perimeter method
        Newton, 1666            2       series related to arcsin x
        Sharp, 1705             1.585   Gregory(sqrt(3))
        Machin, 1706            3.589   Gregory(5, 239)
        Stoermer, 1896          3.168   Gregory(8, 57, 239)

Here Gregory(a1,a2,...,an) is the formula pi/4 = a linear integer
combination of the arctan(1/ai)'s, e.g.  pi/4 = 6 arctan(1/8) + 2
arctan(1/57) + 4 arctan(1/239).  Note that both Sharp's and Stoermer's
"improvements" weren't, showing what happens when you don't do an
explicit complexity analysis.  Both Sharp's and Machin's methods were
used in early computer evaluations of pi.  For some reason Stoermer's
method was later chosen by Shanks and Wrench over Machin's, perhaps due
to inappropriate focus on the dominant arctan: arctan(1/8) must have
seemed more promising than arctan(1/5).

What I want to do here is tie together the first three of these
methods, and connect them up with the quadratrix of Hippias, in a
strikingly simple way.  Oddly enough none of the extant treatments of
pi seem to be aware of this unification, leading one to speculate that
any previous discoverers must have decided it was too simple to bother
mentioning.

The essence of all three methods resides in the following succinctly
stated process step for turning one isosceles triangle into another.

    Make one half of an isosceles triangle isosceles by shortening
    its hypotenuse.

In more detail, in the following figure the initial isosceles triangle
is ROZ, with |OR| = |OZ|.  Take Z' to be the midpoint of RZ, and choose
R' on OR (the hypotenuse of ORZ') so that |OR'| = |OZ'|, yielding a new
isosceles triangle R'OZ'.

          O
             /|\
            / | \
           /  |  \
          /   |   \       Isosceles triangle ROZ becomes isosceles
         /    |    \      triangle R'OZ'.
        R'    |     \
       R------Z'-----Z

Now here's the beautiful part.  This step clearly halves the altitude
from Z to OR, call it y, and also the angle ROZ, call it t, making the
ratio y/t an invariant of the process.  But if we now iterate the
process, in the limit y/t must equal the limiting length of OR, call it
r.  Hence from the initial altitude y0 and the limit r of |OR| we may
recover the initial angle t0 as y0/r.  In particular if we start with
t0 = pi/2, and y0 = 1, then in the limit r = 2/pi.

This is the heart of Archimedes' method of computing pi.

Before getting into any more detail let's rotate the diagram 3/8 of a
turn counterclockwise and coordinatize with X and Y axes thus.

        Y-----Z
            |    /|\
            |  r/ | \
            Y'-/--+--Z'
            | /\  |  |\
            |/ t| |  | \
    --------O-----X--X'LR
        |\_x_/  / e' \___ |z|
        | \_x'_/

So the essence of the process becomes a pair of points X,R moving
respectively right and left, without ever passing.  Furthermore the
line ZR gets steeper at each step, and the length |ZX| halves at each
step, whence the triangle ZXR shrinks to a point, call this the
limiting point L.  This makes it easy to see that the method
converges.

If we regard the point Z as the complex number z then our basic process
step can be described even more succinctly as z := (z + |z|)/2, call
this ALGORITHM A.

Both rectangular coordinates, z = x+iy and polar, z = r exp(it), will
come in handy.  That is, x=Re(z), y=Im(z), r=|z|, t=atan(y/x) (see
diagram).  Primes x',y', etc. denote the corresponding values after the
step.

Precision.  Set e' = r-x' = |X'R|, the *error*, in the sense of a
(weak) upper bound on |X'L|, the remaining distance X has to travel
after the current step taking it to X'.  Since OZ'R is a right-angled
triangle with altitude y', we have x'e' = y'^2, so at the previous
stage xe = y^2.  It follows that as x converges to its limit and y
continues to halve, e must be reduced by a factor of 4 at each step.
Hence as it approaches the limit the method contributes two bits of
precision to the answer at each iteration.

Cost.  Each iteration can be implemented with real arithmetic via
ALGORITHM B:  x := (x + sqrt(x^2+yy))/2, yy := yy/4; where z = x+iy and
yy is a variable doing the job of y^2.  This costs one multiplication
and one square root, i.e. two "major" operations (we will count
division by 4 as an inexpensive shift).  Hence the method obtains one
bit of precision per major operation.

Initializing z to i (t = 90 degrees, r = 1) corresponds to starting
Archimedes' method with a square.  Archimedes actually started with a
hexagon, for which the appropriate initialization is "1 o'clock" (t =
60 degrees, r = 1).

The points visited by z lie on the *quadratrix*, a curve generally
attributed to Hippias, c. 525 BC, and studied in detail two centuries
later by Deinostratos in its role as a device for squaring the circle.
The quadratrix is given by r sin(t) = at  where a is a suitable
constant, e.g. 2/pi when z is initially i.  That the process tracks the
quadratrix follows very easily from r sin(t) = y, which along with t is
halved at each step.  (Neat, eh?)

In Archimedes' construction two parallel series of progressively finer
polygons converge on a circle from both sides.  We may obtain this
effect by modifying the method to renormalize z at each step, scaling
it by an appropriate real so as to keep |z| constant.  This causes z to
now follow a circle instead of a quadratrix.   This entails rewriting
the basic step as ALGORITHM C: z := (z + r)/2; z := z*r/|z|, where r is
now constant, namely the initial value of |z|.  R no longer moves
inward since the scaling keeps resetting it, whence the common limit L
of X and Z must be R.

The points Z and R in this modified construction correspond to two
consecutive vertices of the inside polygon in Archimedes'
construction.  Scaling those two points by the above scaling factor
(that put Z' back on the circle) yields two consecutive vertices of the
outside polygon, with Z' (as scaled) at their midpoint.

We now kill two birds with the one stone: we give Archimedes' analysis
and Vieta's formula using the same reasoning.

Since R goes nowhere it follows that the cumulative effect of all these
scalings, an infinite product, is to stretch the old OL out to the new
(fixed) OR.  If z is initially i then this stretch must be pi/2.
Vieta's formula expresses 2/pi as sqrt(1/2) * sqrt((sqrt(1/2)+1)/2) *
... where the factor following factor x is sqrt((x+1)/2).  We derive
this by showing that these factors are the reciprocals of the scale
factors at each step.

Let z" denote (z+r)/2, the value of z halfway through the step from z
to z'.  Let the line OZ"Z' meet the vertical line x=r (the tangent to
the circle at R) at Q (so the line is OZ"Z'Q in that order).  Q is
called the *inverse* of Z" in the circle and satisfies |OZ"|*|OQ| =
|OZ'|^2 (see e.g. Coxeter, "Intro. to Geometry", section 6.1 on
inversion, or prove it from scratch using the similarity of the
triangles OZ"R and ORQ, and |OR|=|OZ'|).  That is, Z' is at the
geometric mean of Z" and Q; algebraically, z' = sqrt(z"*q).  By
proportionality this relationship must also hold of their x
coordinates, that is, x' = sqrt(((x+1)/2)*r), which when r=1 becomes
sqrt((x+1)/2).  This can be seen to also be the reciprocal of the
desired scaling factor, giving us Vieta's formula for 2/pi as
promised.  Note that Vieta's factors are exactly the values of x we
encounter as we run Algorithm C.  This suggests
ALGORITHM D: x := sqrt((x+1)/2), y := y/(2*x), as the real-valued
version of Algorithm C in the case r=1.

But this same picture can now be used to understand Archimedes' inner
and outer polygons.  Each factor of Vieta's product shrinks the outer
bounding polygon down to the inner one.  The first factor, sqrt(1/2),
shrinks the square, the next the octagon, etc.  (The factors are larger
(closer to 1) if we start from a hexagon instead of a square, with
their product yielding 3/pi instead of 2/pi.)

Archimedes used the perimeter of each polygon to bound the
circumference above and below.  Since OZ"R is congruent to OX'Z', y
gives half the length of one side of the inner polygon, whence y/x (for
r=1) gives the corresponding distance for the outer one.  These values
can be extracted from Algorithm D and multiplied by twice the number of
sides to give the perimeter.  But this multiplication simply undoes the
division of y by 2 at each step of D.  Hence Archimedes' method has for
its sequence of lower bounds on pi (upper bounds on 2/pi) initial
segments of Vieta's product.  For its upper bounds on pi (lower bounds
on 2/pi) it has the same thing but with the sqrt omitted from the last
factor of that segment, since dividing y by x has the effect of
supplying a second copy of the last factor.  (So this also tells us how
to use Vieta's formula directly to bound pi from above and below using
any given initial segment: the segment itself gives the lower bound
while dropping its last sqrt gives the upper.)

Descartes' method.  The term "half-angle" reminds one of the half-angle
formulas for tan(2t) = 2 tan(t)/(1-tan^2(t)), call this H =
2h/(1-h^2).  This inverts as h = (sqrt(1+H^2)-1)/H.  Initializing h to
1, corresponding to t = pi/4, and iterating h := (sqrt(1+h^2)-1)/h  for
k steps, yields t = pi/2^(k+2).  Since tan(t)/t tends to 1 as t tends
to 0 we may approximate pi by h 2^(k+2) after k steps of this process.
This is essentially Descartes' method of computing pi.

But tan(t) is simply y/x, the value we used in conjunction with
Algorithm D for half a side of the outer polygon.  Hence Descartes'
method yields exactly the same values as Archimedes' method for the
outer polygon (so Descartes produces only upper bounds on pi).

(This account of Descartes' isn't as nice as I'd intended; I'd hoped to
equal the standards of the preceding methods by omitting all
trigonometric identities.  I would be very grateful if someone would do
Descartes "right".)

--
Vaughan Pratt   [e-mail]     exp(n pi sqrt(163)/3) for n=1..30