Beeler, M., Gosper, R.W., and Schroeppel, R. HAKMEM. MIT AI Memo 239, Feb. 29, 1972. Retyped and converted to html ('Web browser format) by Henry Baker, April, 1995.

PI

Previous Up Next

ITEM 136: GAUSSIAN INTEGERS (For use by next item.)

Reference: Hardy and Wright, Theory of Numbers. The Gaussian integers are x+iy where x and y are integers. Unique factorization holds, except for powers of i, and the Gaussian primes are
  1. a+bi if a^2+b^2 is prime and
  2. integer primes that = 3 mod 4.
If N(x+iy) = x^2+y^2, then N(u v) = N(u) N(v). If p is prime and not= 3 mod 4, then p = a^2+b^2 has exactly one solution. If n = 3 mod 4, then n = a^2+b^2 has no solution.

To factor x+iy into Gaussian primes, first factor N(x+iy).

  1. If 2 divides N(x+iy), then 1+i and 1-i divide x+iy. Either factor may be used since i (i-1) = i+1.
  2. if p = 3 mod 4 divides N(x+iy), then p divides x+iy.
  3. If p = 1 mod 4 divides N(x+iy) and p = a^2+b^2, then a+ib or b+ia = i (a-ib) divides x+iy. If both do, then p divides x+iy.

ITEM 137 (Salamin): GENERATION OF ARCTANGENT FORMULAS FOR PI

n1 atan(y1/x1) + n2 atan(y2/x2) + ... = n1 arg(x1+iy1) + n2 arg(x2+iy2) + ...

If each x+iy is factored and the n's chosen so all prime factors except 1+i cancel out, the right hand side is a multiple K of pi/4. Some care is needed because of the multiple valuedness of arg. Then, if K = 0, we get an arctangent identity, otherwise we get a pi formula. In the special case of atan(1/x), factorization of x+i is needed. Then case (B) above can't occur, and in case (C), a+ib and a-ib can't both divide x+i.

Example: 

  2
 8  + 1 = 13 x 5

  2             2
18  + 1 = 13 x 5

  2             3
57  + 1 = 13 x 5  x 2
From this we get the factorization 
 8+i = (3+2i) (2-i)

                   2
18+i = (3-2i) (2-i)  i

                   3
57+i = (3-2i) (2+i)  (1-i)
Since we only care about the phase, multiplication by a positive real number may be ignored below. 
     a       b       c         a-b-c      -a-2b+3c      c  b
(8+i)  (18+i)  (57+i)  = (3+2i)      (2+i)         (1-i)  i
We require a-b-c = 0 and -a-2b+3c = 0, which has the minimal non-trivial solution a = 5, b = 2, c = 3. Then we have 
     5       2       3        3  2
(8+i)  (18+i)  (57+i)  = (1-i)  i
Taking the phase of both sides, we get

5 atan(1/8) + 2 atan(1/18) + 3 atan(1/57) = pi/4.
pi formulas: 
pi/4 = atan(1/2) + atan(1/3)
pi/4 = 2 atan(1/3) + atan(1/7)
pi/4 = 4 atan(1/5) - atan(1/239)
pi/4 = 2 atan(1/4) + atan(1/7) + 2 atan(1/13)
pi/4 = 3 atan(1/4) + atan(1/13) - atan(1/38)
pi/4 = 4 atan(1/5) - atan(1/70) + atan(1/99)
pi/4 = 5 atan(1/8) + 2 atan(1/18) + 3 atan(1/57)
pi/2 = 7 atan(1/4) - 5 atan(1/32) + 3 atan(1/132) - 4 atan(1/378)
This last angle has been measured against the International Standard Platinum-Iridium Right angle and certified adequate for any purpose of the U.S. Government, when used in conjunction with a conscientiously applied program of oral hygiene and regular professional care. 
pi/4 = 7 atan(1/9) + atan(1/32) - 2 atan(1/132) - 2 atan(1/378)
pi/4 = 7 atan(1/13) + 8 atan(1/32) - 2 atan(1/132) + 5 atan(1/378)
There are many easily found arctangent identities. Some are: 
  atan(1/31) = atan(1/57) + atan(1/68)
             = atan(1/44) + atan(1/105)
  atan(1/50) = atan(1/91) + atan(1/111)
 atan(1/239) = atan(1/70) - atan(1/99)
             = atan(1/408) + atan(1/577)
atan(1/2441) = atan(1/1164) - atan(1/2225)
             = atan(1/4774) + atan(1/4995)
  atan(1/32) = atan(1/38) + atan(1/132) - atan(1/378)
             = 2 atan(1/73) + atan(1/239) - atan(1/2943)
Infinite sets of arctangent identities: 
     1          1            1
atan(-) - atan(---) = atan(------)
     n         n+1          2
                           n +n+1
Let x(0) = 1, y(0) = 0, x(n) = x(n-1) + 2 y(n-1), y(n) = x(n-1) + y(n-1).

x(n)/y(n) are the continued fraction approximants to sqrt(2).

       1             1              1
atan(-----) + atan(-----) = atan(-------)
     y(2n)         x(2n)         x(2n-1)

       1             1              1
atan(-----) - atan(-----) = atan(-------)
     y(2n)         x(2n)         x(2n+1)

ITEM 138 (Gosper):

pi = 28 arctan (3/79) + 20 arctan (29/278)

pi = 48 arctan (3/79) + 20 arctan (1459/22049)

Which isn't too interesting except that it means that 

       48              20
(79+3i)   (22049+1457i)
is a negative real number.

ITEM 139 (Ramanujan):

 4
 -- =
 pi

 inf
 ====     N
 \    (-1)  (1123 + 21460 N) (1 3 5 ... (2N-1)) (1 3 5 ... (4N-1))
  >   ------------------------------------------------------------
 /                           2N+1   N   3
 ====                     882     32  N!
 N=0
This series gives about 6 decimal places accuracy per term. 
     1
 ---------- =
 sqrt(8) pi

 inf
 ====
 \     (1103 + 26390 N) (1 3 5 ... (2N-1)) (1 3 5 ... (4N-1))
  >    ------------------------------------------------------
 /                        4N+2   N   3
 ====                   99     32  N!
 N=0
This series gives about 8 decimal places accuracy per term. For other pi series, see Ramanujan's paper "Modular Equations and Approximations to Pi" in Quarterly Journal of Pure and Applied Mathematics, vol. 45, page 350 (1914). For more goodies, see "Collected Papers of Srinivasa Ramanujan", Cambridge U. Press (1927).

ITEM 140:

Counting the initial 3 as the zeroth, the 431st denominator in the regular continued fraction for pi is 20776. (Choong, Daykin & Rathbone, Math. of Computation 25 (1971) p. 387).

(Gosper) In the first 26491 terms of pi, the only other 5 digit terms are the 15543rd = 19055 and the 23398th = 19308. (Computed from 35570 terms of the (nonregular) fraction for 4 arctan 1.)

ITEM 141:

The fraction part of 10^760 pi begins: .49999998...

ITEM 142 (Salamin):

Some super-fast convergents to pi if one already has a super-fast computation of trig functions.

X approx pi:

X approx pi/2:

ITEM 143 (Salamin):

Computation of elliptic integrals, log, and pi.

REFERENCES

1. ELLIPTIC INTEGRALS

Define elliptic integrals: 
         1
        /
        [              1
 K(m) = I  ------------------------- dt
        ]             2          2
        /  sqrt((1 - t ) (1 - m t ))
         0

 K'(m) = K(1 - m).
If A(0) and B(0) are given, and then define

AGM(A(0),B(0)) = lim A(n) = lim B(n)

This is called the arithmetic-geometric mean.

Quadratic convergence rate: 

                             2
                  (A(n)-B(n))
A(n+1) - B(n+1) = -----------
                    8A(n+1)
It is known that 
    2              pi
K'(x ) AGM(1, x) = --    [see A&S]
                   2
This gives a super fast method of computing elliptic integrals. It is easy to compute AGM(1, x) for x in the complex plane cut from zero to infinity along the negative real axis. So K'(m) can be computed for -2pi < arg(m) < 2pi, which covers the complex m-plane twice. Handling the phase when taking square roots will permit exploration of more of the Riemann surface.

2. LOGARITHMS

For small m,

K(m) = (pi/2) (1 + m/4 + O(m^2))

e^-pi(K'(m)/K(m)) = (m/16) (1 + m/2 + O(m^2))

Solve for K'(m) and let m = 16/x^2.

K'(16/x^2) = log x + (4/x^2) (log x - 1) + O(log x/x^4).

For x sufficiently large,

log x = K'(16/x^2) = pi/(2 AGM(1, 4/x)).

Requiring a given number of bits accuracy in log x is equivalent to requiring

|K'(16/x^2) - log x)/log x| < epsilon

this becomes

|(4/x^2) (1 - 1/log x)| < |4/x^2| < epsilon

|x| > 2/sqrt(epsilon).

x can be complex. If |x| is not too close to 1, x can be brought into range by reciprocating or repeated squaring.

3. PI

let x = e^n, then

pi = 2 n AGM(1, 4 e^-n).

Suppose epsilon = 10 to the minus a billion.

Than the above equation for pi is valid when n > 1.15 billion.

e^-n is calculated by starting with 1/e and squaring k times.

Thus n = 2^k. 2^30 = 1.07 billion and 2^31 = 2.15 billion, so k = 30 gives 0.93 billion places accuracy and k = 31 gives 1.86 billion places.

ITEM 144 (Schroeppel):

In the above, instead of x = e^n, use x = 2^n and x = e*2^n. Then simultaneous equations can be solve to give both pi and log 2. This avoids having to square e, but requires two AGM's, and therefore takes longer.

Previous Up Next