/*                                                                     
            		    SLOWPI
  	                 PI COMPUTATION                                  
                                                                      
                      Gary M. Blaine 1990                                  

	      Ported to the Amiga by David Hamm
			(March 1991)       
			                                        
    This program is designed to calculate the transcendental          
 constant 'PI' to some arbitrary number of decimal places.  This      
 calculation relies on the fact that 'PI' is exactly equal to:        
                                                                      
       24*arctan(1/8) + 8*arctan(1/57) + 4*arctan(1/239)              
                                                                      
    The program uses a technique whereby extreme precision is         
 obtained through the use of arrays representing "long numbers".      
 Each element of such an array is an integer and contains some        
 number of digits, depending on the word size of the machine,         
 which represent a portion of the "long number".  For example:        

       the number 12.34567890987654321                                

       could be represented in an array as:                           
          array(0) = 12                                               
          array(1) = 345                                              
          array(2) = 678                                              
          array(3) = 909                                              
          array(4) = 876                                              
          array(5) = 543                                              
          array(6) = 210                                              
       if the maximum number of digits per word was equal to 3.       
 Note that the user must interpret where to place the decimal         
 point by some previous convention.                                   

    Since the system arctan function does not deal with "long         
 numbers", this program does.  The subroutines do the actual          
 arithmetic operations of multiplication, division, addition, and     
 subtraction.  Note that it is only necessary for our purposes        
 to do multiplication and division of a "long number" by a            
 "standard number".  The subroutines are restricted in that respect.  

    The arctan function is computed as:                               
                                                                      
       arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 - ...            
                                                                      
 Recall that for an alternating series, the error is less than        
 the absolute value of the last term not taken.                       

    tmp1, tmp2, arctan, and pi must be DIMensioned to the maximum     
 number of words needed to calculate 'PI' to the desired precision    
 using DIGPERWRD digits per word.                                     

 Thanks to Brian Redman for versions of this program in "C" and       
    in "FORTRAN".                                                     
 Brian credits Steve Wampler for the expression calculating the       
    the approximate number of terms to use in computing the arctan    
    function to provide the necessary precision.                      
 All of which required some work to actually get working right.       

 This author has been quite liberal in providing guard digits and     
 extra terms in the computation.                                      

 No attempt has been made to "round" the result.  If you desire       
 to do so, redirect the output to a file and compute to one more      
 place than necessary and manually round with your favorite editor.   

 If compiling using Microsoft C ver. 5.10, use /F 6000 option to       
 increase the stack size.  Smaller values may also work.              

 This version accepts a single command line argument, the number of   
 digits to compute.  If too many or too few arguments are given, the  
 user will be prompted for an input.
*/

#include <math.h>
#include <dos.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define  TRUE        1
#define  FALSE       0
#define  DIGPERWRD   3
#define  MAXNUM   1000
#define  DIM      1500

void set(int number[], int value);
void add(int opr1[], int opr2[], int result[]);
void sub(int opr1[], int opr2[], int result[]);
void mult(int opr1[], int opr2, int result[]);
void divide(int numer[], int denom, int result[]);
void fixit(int number[]);

int nw;
int tmp1[DIM], tmp2[DIM], arctan[DIM], pi[DIM];

void  main(int argc, char *argv[])
{
	char     nu;
	int      np, i, j, nt;
	int      odd;
	unsigned int tstart[2];
	unsigned int tfinish[2];     /* start & finish time value */
	long timeval;           /* time value */
	
	static const int  invx[] = {  8, 57, 239 };
	static const int  coef[] = { 24,  8,   4 };

  printf("SLOWPI v.005\n by Gary M. Blaine and David F. Hamm - Mar 1991\n");

   if(argc != 2)
      {
      fprintf(stderr, "USAGE: SLOWPI {>ram:fileNNN} [number of places]\n ");
      exit(1);
      }
      
   else
      nu = *argv[1];


   if( (nu>'a') && (nu<'Z'))
      {
      fprintf(stderr, "USAGE: SLOWPI {>ram:fileNNN} [number of places]\n");
      exit(1);
      }
      
   else
      np = atoi(argv[1]);

   nw = np / DIGPERWRD + 8;     /* plus some arbitrary fudge */

   if( (np<1) || (nw>DIM) )
      {
      fprintf(stderr, "%d is not a good number\n", np);
      exit(1);
      }

	timeval=timer(tstart);


   set(pi, 0);          /* start with zero */
   for(i=0; i<3; i++)   /* for each arctan term */
      {
      /* compute # terms */
      nt = (int)(.5/log10((double)invx[i])*np)+10;  /* plus fudge */

      /* do first term manually */
      set(tmp1, 1);
      divide(tmp1, invx[i], tmp1);
      set(arctan, 0);
      add(arctan, tmp1, arctan);
      
      odd = TRUE;

      /* compute remaining terms */
      for(j=2; j<=nt; j++)
         {
         odd = !odd;
         divide(tmp1, invx[i], tmp1);
         divide(tmp1, invx[i], tmp1);
         divide(tmp1, 2*j-1,   tmp2);
         if(odd)
            add(arctan, tmp2, arctan);
         else
            sub(arctan, tmp2, arctan);
         }

      /* multiply by coefficient and add */
      mult(arctan, coef[i], arctan);
      add(pi, arctan, pi);
      }

	timeval=timer(tfinish);

   printf("\nTime to compute PI to %d places is %d seconds.\n",
           np, (tfinish[0]-tstart[0]));

   printf("\n\n\n");
   printf("                       THE VALUE OF PI IS APPROXIMATELY:\n");
   printf("\n\n                                     3. +\n\n");

   for(i=0; i<np; i++)
      {
      pi[0] = 0;           /* zero out integer part            */
      mult(pi, 10, pi);    /* multiply by 10 to get next digit */

      if(i % 5 == 0)
         printf(" ");

      if(i % 50 == 0)
         printf("\n         ");

      printf("%c", (char)(pi[0] + '0'));
      }
   
   printf("\n\n\n     ... with accuracy to %d digits\n", np);
   exit(0);
}


/* set a "long number" to an integer value (typically, 0 or 1) */

void set(int number[], int value)
{
   int   i;

   number[0] = value;

   for(i=1; i<nw; i++)
      number[i] = 0;
}


/* add two "long numbers" */

void add(int opr1[], int opr2[], int result[])
{
   int   i;

   for(i=0; i<nw; i++)
      result[i] = opr1[i] + opr2[i];
   
   fixit(result);  /* propagate carries, sort of.  "decimal adjust" */
}


/* subtract two "long numbers" */
/* the original routine had an extra, unnecessary term in the computation */

void sub(int opr1[], int opr2[], int result[])
{
   int   i;

   for(i=0; i<nw; i++)
      result[i] = opr1[i] - opr2[i];
   
   for(i=nw-1; i>=0; i--)
      if(result[i]<0)
         {
         result[i] = result[i] + MAXNUM;   /* propagate borrows */
         result[i-1] = result[i-1] - 1;
         }
}


/* multiply "long number" by small integer...small->try not to overflow */

void mult(int opr1[], int opr2, int result[])
{
   int   i;

   for(i=nw-1; i>=0; i--)
      result[i] = opr1[i] * opr2;
   
   fixit(result);  /* propagate "carries" */
}


/* divide "long number" by integer... avoid overflowing remainder         */
/* casts to long may look like over-kill, but at least some of it is      */
/* absolutely necessary                                                   */

void divide(int numer[], int denom, int result[])
{
   int   temp[DIM], i;
   long  remain;

   remain = 0;
   for(i=0; i<nw; i++)
      {
      temp[i] = (int)(((long)numer[i] + remain)/(long)denom);
      remain=((remain+(long)numer[i])-(long)temp[i]*(long)denom)*(long)MAXNUM;
      }

   for(i=0; i<nw; i++)
      result[i] = temp[i];
}


/* fixup of "long number" so each word is less than MAXNUM */

void fixit(int number[])
{
   int   i;

   for(i=nw-1; i>0; i--)
      if(number[i]>=MAXNUM)
         {
         number[i-1] = number[i-1] + number[i]/MAXNUM;
         number[i] = number[i] % MAXNUM;
         }
}