greg@garnet.berkeley.edu (Greg Kuperberg) (04/03/91)
In article <1991Apr3.014832.15021@linus.mitre.org> bs@faron.mitre.org (Robert D. Silverman) writes: >In article <AVENGER.91Apr2185942@wpi.WPI.EDU> avenger@wpi.WPI.EDU (Samuel Joseph Pullara) writes: >>Could someone please send me a program that will calculate PI to any >>arbitrary number of digits? Thanks in advance... >Let me add my opinion. >Please do NOT comply with the above request. That's a great idea! Here is a simple C program to compute n digits of Pi in time O(n^2). The algorithm is based on the identity atan(1) = 2*atan(1/3)+atan(1/7) and the Taylor series for arctangent. The typical term of arctangent, say 1/2^n/n for n odd, is computed in decimal by base conversion from its representation in a "base" for which the first digit is base n and the rest are base 2 (actually base 4 is more convenient). All terms are computed simultaneously to save the algorithm from cubic running time. Enjoy the program! ---------------------Cut here------------------------------------- /* A program to compute Pi to an arbitrary number of digits. */ /* cc -o pi pi.c */ #define NUM 300 /* Number of BASE digits computed */ #define EXTRA 3 /* Extra digits not printed out */ #define NUMI 1450 /* Space for internal digits. (>NUM*9/8*LOGBASE) */ #define SAFETY 20 /* Extra expansion so that digits aren't dropped */ #define BASE 10000 /* Base we're working in */ #define LOGBASE 4 /* Its logarithm */ #define DISPLAY "%04d" /* For output */ long int d[NUM+1]; void katan1byn(k,n) /* atan(x) = sum of (-1)^i*1/n^(2i+1)/(2i+1) */ int k,n; { long int i,j,j9,n2,c,t,h[NUMI+SAFETY],in[NUMI+SAFETY]; n2 = n*n; in[0] = k*n; for(i=1;i<NUMI+SAFETY;i++) in[i] = h[i] = 0; /* Converts .0...0k from base n...n(2*i+1) to decimal */ for(j=0;j<NUM;j++) { c = t = 0; j9 = (j*9*LOGBASE)/8+SAFETY; for(i=0;i<j9;i++) { in[i] = in[i]*BASE+c; c = in[i]/n2; in[i] %= n2; h[i] = h[i]*BASE+c; t += i&1?-(h[i]/(i+i+1)):h[i]/(i+i+1); h[i] %= i+i+1; } for(i=j+1;;i--) { if(!t) break; d[i] = d[i]+t; c = d[i]%BASE; if(c<0)c+=BASE; t = (d[i]-c)/BASE; d[i] = c; } } } main() { long int i; for(i=0;i<NUM+1;i++) d[i] = 0; katan1byn(8,3); katan1byn(4,7); /* pi = 8*atan(1/3)+4*atan(1/7) */ printf("%1d.",d[0]); for(i=1;i<NUM+1-EXTRA;i++) printf(DISPLAY,d[i]); printf("\n"); return 0; } ---- Greg Kuperberg greg@math.berkeley.edu
jroth@allvax.enet.dec.com (Jim Roth) (04/04/91)
In article <1991Apr3.070614.11682@agate.berkeley.edu>, greg@garnet.berkeley.edu (Greg Kuperberg) writes... >In article <1991Apr3.014832.15021@linus.mitre.org> bs@faron.mitre.org (Robert D. Silverman) writes: >>In article <AVENGER.91Apr2185942@wpi.WPI.EDU> avenger@wpi.WPI.EDU (Samuel Joseph Pullara) writes: >>>Could someone please send me a program that will calculate PI to any >>>arbitrary number of digits? Thanks in advance... >>Let me add my opinion. >>Please do NOT comply with the above request. > >That's a great idea! Here is a simple C program to compute n digits of >Pi in time O(n^2). Here's an even shorter program :-) Stan and I once sat down and wrote a one line TECO macro that did this calculation! /* 1000 digits of PI */ /* 'Spigot' algorithm origionally due to Stanly Rabinowitz */ #include stdio #include math 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"); } } >Enjoy the program! - Jim