Pi
In high school, I wrote a program to generate pi (3.14159...) to 75,000 digits.
I was rather interested in the subject and tried several different methods
of generating the number.
Here is the output of the program separated into lines
(here without lines).
Some of the formulas I used were the following:
π = 2 arcsin 1
π/2 = 1 + 1/(2*3) + 3/(2*4*5) + (3*5)/(2*4*6*7) + (3*5*7)/(2*4*6*8*9)
+ ...
π = 6 arcsin ½
π/6 = 1/2 + 1/(2*3*2^3) + 3/(2*4*5*2^5) + (3*5)/(2*4*6*7*2^7) + ...
arctan x = Σ x(-x²)ⁿ/(2n+1)
Note: where y = x²/(x²+1),
arctan x = (y/x)(1 + (2/3)y + (2*4/3*5)y^2 + (2*4*6/3*5*7)y^3 + ...)
π = 20 arctan (1/7) + 8 arctan (3/79)
π = 4 arctan 1 (Slow!)
π/4 = 1/1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + 1/13 - 1/15 + ...
(2 decimal places cannot be reached with 300 terms)
π/4 = arctan(1/2) + arctan(1/5) + arctan(1/8)
π = 16 arctan (1/5) - 4 arctan (1/239)
π = 24 arctan (1/8) + 8 arctan (1/57) + 4 arctan (1/239)
π/4 = arctan (1/2) + arctan (1/8) + arctan (1/18) + ...
( ... = arctan (13/91) )
n sin Θ < π < n tan Θ
( Inscribed and circumscribed figures )
⌠¼
π = (3/4)√3 + 24 │ √(x - x²) dx
(Fast!)
⌡0
π = (3/4)√3 + 24 Σ [(-1)*1*3*5...(2i-3)]/[2*4*6*8...(2i)](2i+3)2^(2i+2)
π = (3/4)√3 + 24(1/(3*2^2) - 1/(2)(5*2^4) - 1/(2*4)(7*2^6) - (1*3)/(2*4*6)(9*2^8)
- ...)
pi = (3/4)√3 + 2 ; with first term already
n = 2 ; 2i
a = 24 ; -24(numerator)
b = 2*2^4 ; (denominator)/(2i+3)
┌ pi -= a/(n+3)b ; subtract next term
│ n += 2 ; increment i
│ a *= (n-3) ; next factor of numerator
│ b *= n*(2^2) ; next factor of denominator
└─loop
Alternatively,
pi = (3/4)√3 + 2 ; with first term already
n = 2 ; 2i
t = 24/(2*2^4) ; -24(2i+3)(ratio)
┌ pi -= t/(n+3) ; subtract next term
│ n += 2 ; increment i
│ t *= (n-3)/4n ; factors of next term
└─loop
Accepted:
π²/6 = 1/1² + 1/2² + 1/3² + 1/4² + 1/5²
+ ... (Slow!)
π²/8 = 1/1² + 1/3² + 1/5² + 1/7² + ...
π²/12 = 1/1² - 1/2² + 1/3² - 1/4² + ...
π = 4/(1 + 1²/(2 + 3²/(2 + 5²/(2 + 7²/(2 + 9²/(2
+ 11²/(...))))))) (Slow!)
π = 2/√½*√(½+½√½)*√(½+½√(½+½√½))*...
x = 0 ; cos 90°
pi = 2 ; numerator
┌ x = √(½ + ½x) ; x = cos ½Θ
│ pi /= x ; multiply denominator
└─loop
Unsure:
π²/6 = (2²/(2²-1)) * (3²/(3²-1)) * (5²/(5²-1))
* (7²/(7²-1)) * ...
π = (2/1)(2/1)(2/3)(4/3)(4/5)(6/5)(6/7)(8/7)(8/9)(10/9)(10/11) ...
I had my pi and e generators all ready to send you, and then I remembered
that I
was going through the Route Y mailer, which won't do attachments. I'll
have to
bring them on a disk to school and e-mail them to you from a school computer.
Until then, here's the formulas (formulae) about pi that I had lying around.
If
the math symbols come through all garbled or not there at all, then tell
me and
we'll change the way the formulas are written.
These are the formulas that I gathered as I was trying to generate pi:
Two formulas with the arcsin, evaluated with the power series for arcsin.
pi/2 = arcsin 1 = 1 + 1/(2*3) + 3/(2*4*5) + (3*5)/(2*4*6*7) +
(3*5*7)/(2*4*6*8*9) + ...
pi/6 = arcsin ½ = 1/2 + 1/(2*3*2^3) + 3/(2*4*5*2^5) +
(3*5)/(2*4*6*7*2^7) + ...
Several series involving the arctan, evaluated with the power series,
arctan x = sum x(-x²)^n/(2n+1)
or Euler's faster method, taking y = x²/(x²+1),
arctan x = (y/x)(1 + (2/3)y + (2*4/3*5)y^2 + (2*4*6/3*5*7)y^3 + ...)
pi = 4 arctan 1 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - 1/11 + ......
pi = 20 arctan (1/7) + 8 arctan (3/79)
pi = 4 arctan(1/2) + 4 arctan(1/5) + 4 arctan(1/8)
pi = 16 arctan (1/5) - 4 arctan (1/239)
pi = 24 arctan (1/8) + 8 arctan (1/57) + 4 arctan (1/239)
pi/4 = arctan (1/2) + arctan (1/8) + arctan (1/18) + arctan (13/91)
Next, there were the inscribed and circumscribed figures. Draw a circle,
divide
the center angle n ways, so t=2pi/n, and then 2 sin t is the length of a
side of
an inscribed regular polygon of n sides, and 2 n tan t is the length of a
side
of a similar circumscribed polygon, so we have
n sin t < pi < n tan t
Next draw a circle with radius ½ and center ½, and integrate
it (sqrt(x - x²))
from 0 to ¼. Multiply by 24 and add (3/4)sqrt3. You get
pi. Then you can
estimate the integral with a sum,
pi = (3/4)sqrt3 + 24 sum [(-1)*1*3*5···(2i-3)]/
[2*4*6*8···(2i)](2i+3)2^(2i+2)
pi = (3/4)sqrt3 + 24(1/(3*2^2) - 1/(2)(5*2^4) - 1/(2*4)(7*2^6) -
(1*3)/(2*4*6)(9*2^8) - ...)
Here's an algorithm to use this method, taking the fractions exactly:
pi = (3/4)sqrt3 + 2 ; with first term already
n = 2 ; 2i
a = 24 ; -24(numerator)
b = 2*2^4 ; (denominator)/(2i+3)
(top of loop)
pi = pi - a/(n+3)b ; subtract next term
n = n + 2 ; increment i
a = a(n-3) ; next factor of numerator
b = bn(2^2) ; next factor of denominator
loop
Alternatively, using a floating point for the term,
pi = (3/4)sqrt3 + 2 ; with first term already
n = 2 ; 2i
t = 24/(2*2^4) ; -24(2i+3)(ratio)
(top of loop)
pi = pi - t/(n+3) ; subtract next term
n = n + 2 ; increment i
t = t(n-3)/4n ; factors of next term
loop
Here's a few others. These next few deal intimately with
the Riemann Zeta
Function, which is what my current research project is about. Unfortunately,
as
algorithms, they're pretty slow.
pi²/6 = 1/1² + 1/2² + 1/3² + 1/4² + 1/5² +
...
pi²/8 = 1/1² + 1/3² + 1/5² + 1/7² + ...
pi²/12 = 1/1² - 1/2² + 1/3² - 1/4² + ...
These two are a little odd; each is an infinite product, and also pretty
slow:
pi²/6 = (2²/(2²-1)) * (3²/(3²-1)) * (5²/(5²-1))
* (7²/(7²-1)) * ...
pi = (2/1)(2/1)(2/3)(4/3)(4/5)(6/5)(6/7)(8/7)(8/9)(10/9)(10/11) ...
Here is a continued fraction for pi. I hear that continued fractions
are all
the rage in numerical methods (computer calculation), but this particular
one is
quite slow.
pi = 4/(1 + 1²/(2 + 3²/(2 + 5²/(2 + 7²/(2 + 9²/(2
+ 11²/(...)))))))
Another: 2/pi = cos 45 cos 22.5 cos 11.25 ...
pi=2/sqrt½*sqrt(½+½sqrt½)*sqrt(½+½sqrt(½+½sqrt½))*...
x = 0 ; cos 90
pi = 2 ; numerator
(loop top)
x = sqrt(½ + ½x) ; x = cos
½t
pi /= x ; multiply denominator
loop
Those were the ones I had gathered long ago. Recently, looking for
libraries
for C that process numbers of arbitrary precision, I ran into some more.
Here
are the new things I found: One was a little C program to generate pi.
I also
found another pair of C programs to generate pi and e, respectively.
I'll only
include the one about pi, because there's a bunch more formulas for e.
-------------------------------------------------------------
In the latest revision of Numerical Recipes, there is a Brent-Salamin
AGM algorithm to compute PI, but I don't have it online.
For just 10K digits of PI, a simple program like this will suffice.
- Jim
/* 1000 digits of PI
*/
/* 'Spigot' algorithm origionally due to Stanly Rabinowitz */
#include <stdio.h>
main()
{
int d = 4, r = 10000, n = 251, m = 3.322*n*d;
int i, j, k, q;
static int a[3340];
for (i = 0; i <= m; i++) a[i] = 2;
a[m] = 4;
for (i = 1; i <= n; i++) {
q = 0;
for (k = m; k > 0; k--) {
a[k] = a[k]*r+q;
q = a[k]/(2*k+1);
a[k] -= (2*k+1)*q;
q *= k;
}
a[0] = a[0]*r+q;
q = a[0]/r;
a[0] -= q*r;
printf("%04d%s",q, i & 7 ? " " : "\n");
}
}
int a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}
-------------------------------------------------------------
Said someone else in the same document,
-------------------------------------------------------------
Providing you have an efficient subroutine package for multiprecision
arithmetic, there are a number of quickly convergent iterations which
can be employed to calculate many digits of Pi.
Before the latest approximations to Pi using modular equations were
developed, Taylor series expansions of arctans were the state of the
art. There are a huge number of these, all designed to make the
denominators of the series terms grow as rapidly as possible.
Besides the one above, the following are also popular...
pi = 20*arctan(1/7) + 8*arctan (3/79)
pi/4 = arctan(1/2) + arctan(1/5) + arctan(1/8)
pi/4 = 3*arctan(1/4) + arctan(1/20) + arctan(1/1985)
pi = 24*arctan(1/8) + 8*arctan(1/57) + 4*arctan(1/239)
pi = 48*arctan(1/18) + 32*arctan(1/57) - 20*arctan(1/239)
pi/4 = 22*arctan(1/28)+2*arctan(1/443)-5*arctan(1/1393)-10*arctan(1/11018)
pi/4 = 44*arctan(1/57)+7*arctan(1/239)-12*arctan(1/682)+24*arctan(1/12943)
My favorite is a quartic iteration based on an elliptic function called
the Singular Value Function of the Second Kind. It converges exponentially
to 1/Pi as its argument increases. There is a rich set of "modular
equations" which algebraically relate values of this function for
different values of its argument. Repeating these formulas as iterations
allows us to compute the value of this function for increasingly large
arguments and consequently the value of Pi.
This technique is based on the work of the eccentric Indian
mathematician Srinivasa Ramanujan.
The basic iteration goes like this...
Y[N+1] = (1-(1-Y[N]^4)^0.25)/(1+(1-Y[N]^4)^0.25)
A[N+1] = (1+Y[N+1])^4*A[N]-4^(N+1)*SQRT(R)*Y[N+1]*(1+Y[N+1]+Y[N+1]^2)
This will compute approximately 45 million decimal digits of Pi in
only 12 iterations with R = 4.
This formula works for any positive R, but the initial values of the
iteration, Y[0] and A[0] are dependent on the R chosen. A[0] should
be the Singular Value Function of the Second Kind evaluated at R and
Y[0] should be the square root of the another function, called the
Singular Value Function of the First Kind, again evaluated at R.
For instance, if we choose R = 2, the correct starting values turn out
to be A[0] = SQRT(2)-1 and Y[0] = SQRT(SQRT(2)-1).
Then the iteration looks like this...
N Y[N] A[N]
1/A[N]
0 0.6435942530 0.414213562 2.41421356
1 0.0235239602 0.318310704 3.14158459
2 0.0000000383 0.318309886 3.14159265
As you can see, the convergence to Pi is quite rapid, and each iteration
gives approximately four times the number of decimal places as the
previous one.
Using techniques like these, and an efficient multiprecision package,
computation of Pi to millions of places is well within the computational
power of the typical personal computer.
---------------------------------------------------------------
Francis Barrett, F.R.C. | Thou canst not travel on the path |
The Cybernetics Guild | before thou hast become the Path
|
fb@cyberg.win.net | itself.
|
---------------------------------------------------------------
Finally, I'll end with some evaluations I've done in hexadecimal and binary,
and
then I'll get my pi and e generators to you from another computer.
Hexadecimal approximation:
π ≈ 3.243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA
98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F
84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADF
B7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E285
8EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD
5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA
417918B8DB38EF8E79DC
Binary approximation:
π ≈ 11.00100100001111110110101010001000100001011010001100001
000110100110001001100011001100010100010111000000011011100000
111001101000100101001000000100100111000001000100010100110011
111001100011101000000001000001011101111101010011000111011000
100111001101100100010010100010100101000001000011110011000111
000110100000001001101110111101111100101010001100110110011110
011010011101001000011000110110011000000101011000010100110110
111110010010111110001010000110111010011111110000100110101011
011010110110101010001110000100100010111100100100001011011010
101110110011000100101111001111110110001101111010001001100010
000101110100110100110001101111110110101101011000010111111111
101011100101101101111010000000110101101111110110111101110001
110000110101111111011010110101000100110011111101001011010111
010011111001001000001000101111100010010110001111111100110010
010010010100001100110010100011110110011100100010110110011110
111000010000000000111110010111000101000010110001110111111000
001011001100011011010010010000011011000011100010101011101001
110011010011010010001011000111111101010001111110100100100110
011110101111110000011011001010101110100100011110111001010001
110101101100101100001110001100010111100110101011000100000100
001010101001010111011100111101101010100101001000001110111000
010010110100101100110110101100111000011000011010101001110010
010101011110010011000000001001111000101110100011011000000100
011001010000110000010000101111100001100101001000001011110010
001100010111000110110110011100011101111100011100111100111011
1001
Here is a hexadecimal approximation to pi:
π ≈ 3.243F6A8885A308D313198A2E03707344A4093822299F31D0082EF
A98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C
50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FF
D72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B
3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D
95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392A
F26013C5D1B023286085F0CA417918B8DB38EF8E79DC
Here is a binary approximation to pi:
π ≈ 11.0010010000111111011010101000100010000101101000110000100
0110100110001001100011001100010100010111000000011011100000111
0011010001001010010000001001001110000010001000101001100111110
0110001110100000000100000101110111110101001100011101100010011
1001101100100010010100010100101000001000011110011000111000110
1000000010011011101111011111001010100011001101100111100110100
1110100100001100011011001100000010101100001010011011011111001
0010111110001010000110111010011111110000100110101011011010110
1101010100011100001001000101111001001000010110110101011101100
1100010010111100111111011000110111101000100110001000010111010
0110100110001101111110110101101011000010111111111101011100101
1011011110100000001101011011111101101111011100011100001101011
1111101101011010100010011001111110100101101011101001111100100
1000001000101111100010010110001111111100110010010010010100001
1001100101000111101100111001000101101100111101110000100000000
0011111001011100010100001011000111011111100000101100110001101
1010010010000011011000011100010101011101001110011010011010010
0010110001111111010100011111101001001001100111101011111100000
1101100101010111010010001111011100101000111010110110010110000
1110001100010111100110101011000100000100001010101001010111011
1001111011010101001010010000011101110000100101101001011001101
1010110011100001100001101010100111001001010101111001001100000
0001001111000101110100011011000000100011001010000110000010000
1011111000011001010010000010111100100011000101110001101101100
111000111011111000111001111001110111001