/*
 *  This program calculates the value of PI to many decimal places.
 *  It is patterned after Bob Bishop's "Apple Pi" program which was
 *  written in BASIC for the Apple II computer.  This version uses
 *  a fixed memory size for the calculation arrays defined by
 *  MEMSIZE.  Here MEMSIZE is set to 800,000, which means the maximum
 *  number of decimal digits that can be calculated will be 
 *  approximately 1.6 million.  To enable the calculation of more
 *  decimal digits of Pi, just increase the value of MEMSIZE.
 *
 */

#include  <stdio.h>
#include  <stdlib.h>
#include  <string.h>
#include  <time.h>

#define  TEN         100
#define  MEMSIZE     800000

unsigned char  power[MEMSIZE];
unsigned char  term[MEMSIZE];
unsigned char  result[MEMSIZE];

int  constant[] = {25,239};

int  coef, sign, size, zero;

void  add();
void  calc_pi();
void  copy_array();
void  divide();
void  init_calc();
void  print_result();
void  subtract();

/*------------------------------------------------------------------------*/

int main(argc,argv)
int   argc;
char  *argv[];
{
   char  chsize[24];
   int   timing;

   if(argc != 2)
     {
       fprintf(stderr,"--> Usage:  pi  <number>\n");
       exit(-3);
     }

   strcpy(chsize,argv[1]);

   size = atoi(chsize);

   size = (size + 1) / 2;

   size += 5;

   timing = clock();

   calc_pi();

   print_result();

   timing = clock() - timing;

   fprintf(stderr,"Time is %f seconds.\n",(double) timing/CLOCKS_PER_SEC);
}

/*------------------------------------------------------------------------*/

void calc_pi()
{
   int  i = 0, pass;

   for(pass = 0; pass < 2; pass++)
     {
       init_calc(size,pass);

       do
         {
           copy_array(size);

           divide(term,coef);

           if(sign > 0)
               add();
           else if(sign < 0)
               subtract();

           coef += 2;
           sign = -sign;

           divide(power,constant[pass]);

           /* printf("interation %5d, zero = %d\n",i++,zero); */
           /* printf("-- coef = %d, term[0] = %d\n",coef,*term); */

           if(pass == 1)
               divide(power,constant[pass]);
         }
       while(zero);

       /* printf("Next pass...\n"); */
     }  /* end of for(i) */
}

/*------------------------------------------------------------------------*/

void divide(digits,divisor)
unsigned char  digits[];
int  divisor;
{
   unsigned char  *cptr;
   int  digit = 0;
   int  quotient, residue;

   /* printf("divide:  array = %p, divisor = %d\n",digits,divisor); */
   /* printf("divide:  number starts with %d\n",*digits); */

   zero = 0;

   for(cptr = digits; cptr < (digits+size); cptr++)
     {
       digit += *cptr;
       quotient = digit / divisor;
       residue  = digit % divisor;
       zero |= quotient + residue;
       *cptr = quotient;
       digit = TEN * residue;
     }
}

/*------------------------------------------------------------------------*/

void  add()
{
   int  carry, i, sum;

   carry = 0;

   for(i = size - 1; i > -1; i--)
     {
       sum = result[i] + term[i] + carry;
       
       if(sum >= TEN) 
         {
           sum -= TEN;
           carry = 1;
         }
       else
         {
           carry = 0;
         }

      result[i] = sum;
    }  /* end for(i) */
}

/*------------------------------------------------------------------------*/

void  subtract()
{
   int  difference, i, loan;

   loan = 0;

   for(i = size - 1; i > -1; i--)
     {
       difference = result[i] - term[i] - loan;
       
       if(difference < 0) 
         {
           difference += TEN;
           loan = 1;
         }
       else
         {
           loan = 0;
         }

       result[i] = difference;
    }  /* end of for(i) */
}

/*------------------------------------------------------------------------*/

void  init_calc(size,pass)
int  size, pass;
{
   int  divisor, i;

   pass++;

   memset(power,0,size);
   memset(term,0,size);

   if(pass == 1)
       memset(result,0,size);

   *power = 16 / (pass * pass);

   switch(pass)
     {
       case 1:
           divisor = 5;
           sign = 1;
           break;
       case 2:
           divisor = 239;
           sign = -1;
           break;
     }
   
   divide(power,divisor);
   
   coef = 1;

   /* printf("power = %d\n",*power); */
   /* printf("sign  = %d\n",sign); */
   /* printf("pass  = %d\n",pass); */
}

/*------------------------------------------------------------------------*/

void  copy_array(size)
int  size;
{
   memcpy(term,power,size);

   /* fprintf(stderr,"\a\a\a\a\a"); */
}

/*------------------------------------------------------------------------*/

void print_result()
{
   int            cnt;
   unsigned char  *cptr;

   size -= 5;

   printf("The value of pi to %d decimal places is:\n",(TEN/100+1)*size);

   printf("     %d.\n",*result);

   cnt = 0;

   for(cptr = result+1; cptr < result+size; cptr++)
     {
       cnt++;
       printf("%02d",*cptr);

       if(!(cnt % 30))
           printf("\n");
       else if(!(cnt % 3))
           printf("  ");
     }

   printf("\n");
}

/*------------------------------------------------------------------------*/