/*----------------------------------------------------------------------
program to calculate pi, version 5.1f

Almost exactly like version 4.8, except it uses floating point
and works with 9 digit chunks instead of 4 digit chunks. Likewise,
using floating point allows multiplication by reciprocals instead
of (much slower) division. For 5000 digits it needs 6.92 seconds on
my computer; however, this program runs over 30% faster than the
integer version on a Pentium, and seven or eight times faster on a RISC
machine like a Sun. Finally, 53 bit words mean this program is capable
of over 6.3 million digits of pi. 

Version 5.1f is a bug fix of version 5.0f, that simplifies some of the
code (for a 1% speedup over 5.0f) and also fixes some problems with
carries and borrows that put wrong digits all over the computed results.
Up to 32000 digits I'm confident it gives correct results.

Questions, comments, and (especially) speedups should be directed to
jasonp@isr.umd.edu
-----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

void printout( void );
void atan5( void );
void atan239( void );

long *sum, *term, firstword, lastword, words, denom;
double recip, TEN9=1000000000;

int main(int argc, char *argv[]) {

    double dbldiv;
    clock_t end,start;
    long digits = 0, temp, remainder, i, quotient;

    if(argc != 2){
        printf("\nusage: pi50f NumberOfDigits > OutputFile\n\n");
        exit( EXIT_FAILURE );
    }
    else{
        digits = atol( argv[1] );
    }

                                              /* Allocate array space */
   words = digits/9 + 2;
   sum = (long *)calloc( words+2, sizeof(long) );
   term = (long *)calloc( words+2, sizeof(long) );
   if (term == NULL || sum == NULL ) {
       printf("Memory allocation failed. Try fewer digits.\n");
       exit( EXIT_FAILURE );
   }

                                          /* ----- 16*atan(1/5) -------*/ 
   start = clock();                             
   firstword = 0; lastword = 1; denom = 3;
   sum[0] = 3; term[0] = 3; recip = .04;      
   sum[1] = 2e8; term[1]=2e8;

   while( firstword<words ) {
      atan5();
   }


                                        /* ----- -4*atan(1/239) ------- */
   remainder = 4; recip = 1.0 / 239;

   for( i=1; i<=words; i++ ){                    /*crunch out first term*/
      temp = sum[i];
      dbldiv = (double)remainder * TEN9;
      quotient = dbldiv * recip;
      remainder = dbldiv - 239.0 * quotient + .5;
      term[i] = quotient;
      temp -= quotient;
      if ( temp<0 ) { temp += TEN9; --sum[i-1]; }
      sum[i] = temp;
   }

   firstword=1; denom=3; recip = 1.0/57121;

   while(firstword<words){
      atan239();
   }
                                         
   end = clock();                                      /* -----Finish up */
   printout();
   printf("\nComputation time = %6.2f seconds\n",
         (float)(end-start)/(float)CLOCKS_PER_SEC );

   free(sum); free(term); exit(0);
}
/*------------------------------------------------------------------------*/
void atan5( void ) {

   long remainder1=0, remainder2=0, remainder3=0, remainder4=0,
        i, quotient, denom2=denom+2;
   double dbldiv, recipdenom = 1.0/denom, recipdenom2 = 1.0/denom2;

   for( i=firstword; i<=lastword+1; i++) {
                                                            /* minus next */
      dbldiv = remainder1 * TEN9 + term[i];
      term[i] = dbldiv * recip;
      remainder1 = dbldiv - 25.0 * term[i] + .5;

      dbldiv = remainder2 * TEN9 + term[i];
      quotient = dbldiv * recipdenom;
      remainder2 = dbldiv - (double)denom * quotient + .5;
      sum[i] -= quotient;
      if ( sum[i]<0 ) {
         sum[i] += TEN9; --sum[i-1];
         if ( sum[i]<0 ) { sum[i] += TEN9; --sum[i-1]; }
      }
                                                             /* plus next */
      dbldiv = remainder3 * TEN9 + term[i];
      term[i] = dbldiv * recip;
      remainder3 = dbldiv - 25.0 * term[i] + .5;

      dbldiv = remainder4 * TEN9 + term[i];
      quotient = dbldiv * recipdenom2;
      remainder4 = dbldiv - (double)denom2 * quotient + .5;
      sum[i] += quotient;
      if ( sum[i]>=TEN9 ) {
         sum[i] -= TEN9; ++sum[i-1];
         if ( sum[i]>=TEN9 ) { sum[i] -= TEN9; ++sum[i-1]; }
      }
   }
                                                    /*finish off the rest */
   if (term[i-1]>0 && lastword<words) lastword++;   /*of terms; divides by*/
                                                    /*25 don't need full  */
   for(; i<=words; i++) {                           /*precision           */
      dbldiv = remainder2 * TEN9;
      quotient = dbldiv * recipdenom;
      remainder2 = dbldiv - (double)denom * quotient + .5;
      sum[i] -= quotient;
      if ( sum[i]<0 ) {
         sum[i] += TEN9; --sum[i-1];
         if ( sum[i]<0 ) { sum[i] += TEN9; --sum[i-1]; }        
      }

      dbldiv = remainder4 * TEN9;
      quotient = dbldiv * recipdenom2;
      remainder4 = dbldiv - (double)denom2 * quotient + .5;
      sum[i] += quotient;
      if ( sum[i]>=TEN9 ) {
         sum[i] -= TEN9; ++sum[i-1];
         if ( sum[i]>=TEN9 ) { sum[i] -= TEN9; ++sum[i-1]; }
      }
   }

   if( term[firstword] == 0 ) firstword++;
   denom+=4;
}

/*------------------------------------------------------------------------*/
void atan239(void) {

   long remainder1=0, remainder2=0, remainder3=0, remainder4=0,
        i, denom2 = denom + 2, quotient;
   double dbldiv, recipdenom = 1.0/denom, recipdenom2 = 1.0/denom2;

   for( i=firstword; i<=words; i++) {                      /* plus next */
      dbldiv = remainder1 * TEN9 + term[i];
      term[i] = dbldiv * recip;
      remainder1 = dbldiv - 57121.0 * term[i] + .5;

      dbldiv = remainder2 * TEN9 + term[i];
      quotient = dbldiv * recipdenom;
      remainder2 = dbldiv - (double)denom * quotient + .5;
      sum[i] += quotient;
      if ( sum[i]>=TEN9 ) {
         sum[i] -= TEN9; ++sum[i-1];
         if ( sum[i]>=TEN9 ) { sum[i] -= TEN9; ++sum[i-1]; }
      }
                                                           /* minus next */
      dbldiv = remainder3 * TEN9 + term[i];
      term[i] = dbldiv * recip;
      remainder3 = dbldiv - 57121.0 * term[i] + .5;

      dbldiv = remainder4 * TEN9 + term[i];
      quotient = dbldiv * recipdenom2;
      remainder4 = dbldiv - (double)denom2 * quotient + .5;
      sum[i] -= quotient;
      if ( sum[i]<0 ) {
         sum[i] += TEN9; --sum[i-1];
         if ( sum[i]<0 ) { sum[i] += TEN9; --sum[i-1]; }        
      }
   }

   firstword++;
   if( term[firstword]==0 ) firstword++;
   denom += 4;
}
/*----------------------------------------------------------------------*/
void printout(void){
   long i = 1;
   printf("pi=3.\n");

   for(; i<words; i++) { 
      printf("%.9ld ", sum[i]);
      if( i%7 == 0) printf("\n");
   }
}