- a+b
*i*if a^2+b^2 is prime and - integer primes that = 3 mod 4.

To factor x+*i*y into Gaussian primes, first factor N(x+*i*y).

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

If each x+*i*y 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+*i*b and a-*i*b 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 2From 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) iWe 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) iTaking the phase of both sides, we get

5 atan(1/8) + 2 atan(1/18) + 3 atan(1/57) = pi/4.

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+1Let 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)

*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.

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=0This 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=0This series gives about 8 decimal places accuracy per term. For other

(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.)

X approx *pi*:

- X <- X + sin X, epsilon <- epsilon^3/6
- X <- X - tan X, epsilon <- -epsilon^3/3

- X <- X + cos X, epsilon <- epsilon^3/6
- X <- X + cot X, epsilon <- -epsilon^3/3

REFERENCES

- Whittaker & Watson,
*Modern Analysis,*chap. 22. - Abramowitz & Stegun,
*Handbook of Mathematical Functions,*sect. 17.3, 17.6

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

- A(n+1) =
*arithmetic*mean of A(n) and B(n) - B(n+1) =
*geometric*mean of A(n) and B(n)

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] 2This 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 -2

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.

*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.