Xref: utzoo comp.compression:201 sci.math:16348
Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!pacific.mps.ohio-state.edu!linac!att!ucbvax!agate!garnet.berkeley.edu!greg
From: greg@garnet.berkeley.edu (Greg Kuperberg)
Newsgroups: comp.compression,sci.math
Subject: A program for computing Pi
Message-ID: <1991Apr3.070614.11682@agate.berkeley.edu>
Date: 3 Apr 91 07:06:14 GMT
References: <AVENGER.91Apr2185942@wpi.WPI.EDU> <1991Apr3.014832.15021@linus.mitre.org>
Sender: usenet@agate.berkeley.edu (USENET Administrator)
Reply-To: greg@math.berkeley.edu
Followup-To: sci.math
Organization: U.C. Berkeley
Lines: 57

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