Solving Pi

[eeg systems]

	I am posting this as an interesting approach to solving
pi to many decimal places. Originally compiled and run on an Apple
II, it took 40 minutes to solve pi to 1000 places. On a Vax, it took
1.5 seconds! I limited its precision to 1000 places because the only
listing for pi I had available only listed that many places.

	To use, cut below and compile.

	Enjoy.  Bryan Costales ( ...!dual!ptsfa!ski!eeg!bcx )
 
--------------------------< cut here >-------------------------------
/*
** NAME
**	pi -- a program to calculate the value of pi
**
** SYNOPSIS
**        pi {-}digits {t}
**
** DESCRIPTION
**	Pi calculates the value of pi to the number of digits
**	specified. Output is to the standard output, if the
**	digits specifed is negative, output is formatted
**	for attractive display on an 80col printer.
**
**	The optional last argument "t" enables timing of the
**	computational loop. The pitime() function must be supplied
**	by the individual site implimentation.
**
**	This is an adaptation in C of the BASIC program "Apple Pi"
**	by Robert J. Bishop, published in Micro Apple (c) 1981
**	Micro Ink, Inc. That program, when run on an Apple II,
**	took 40 hours to solve pi to 1000 places. Arrays of digits
**	are used to impliment extended precision arithmatic. The
**	program solves pi using the series expansion:
**
**	      infinity                    infinity
**	      ____  16*(-1e(k+1))         ____  4*(-1e(k+1))
**	      \                           \
**	pi =   >    -------------    -     >    ------------
**	      /                           /
**	      ----  (2k-1)*5e(2k-1)       ----  (2k-1)*239e(2k-1)
**	      k = 1                       k = 1
**
** AUTHOR
**	Bryan Costales ( ...!dual!ptsfa!ski!eeg!bcx )
**	pi.c (c) 1985 Bryan Costales
**
** BUGS
**	
*/
 
 
#include <stdio.h>
 
/*
 * The maximum number of digits (plus 5 for rounding).
 */
#define MAX 10005
 
/*
 * Global variables
 */
int 	Sign, 
	Zero, 
	Size, 
 	Pass,	/* two passes */
	Exp,	/* exponent for divide */
	Divide; 
 
/*
 * These character arrays represent the digits for extended
 * precision arithmatic. (These may need to be int on some machines).
 */
#define DIGIT char
DIGIT Power[MAX], Term[MAX], Result[MAX];
 
main(argc, argv)
int  argc;
char *argv[];
{
	static int constant[3]={0,25,239};
	register int i;
	int charcount, printer = 0;
	void exit();	/* for lint */
 
	if(argc == 1)
	{
		(void)fprintf(stderr,"usage: pi {-}digits {t}\n");
		exit(1);
	}
 
	Size = atoi(argv[1]);
	if(Size < 0)
	{
		Size = -Size;
		printer++;
	}
 
	if(Size >= MAX-4 || Size <= 0)
	{
		(void)fprintf(stderr,"'digits' must be 1 thru %d.\n",MAX-5);
		exit(1);
	}
 
	if( (argc==3) && (*argv[2]=='t') )
	{
		(void)fprintf(stderr,"start computation: ");
		pitime();
	}
 
	for(Pass=1; Pass < 3; Pass++)
	{
		init();
 
		do
		{
			copy();
			Divide = Exp;
			div(Term);
			if ( Sign > 0 ) 
				add();
			if ( Sign < 0 ) 
				sub();
			Exp = Exp + 2;
			Sign *= -1;
			Divide = constant[Pass];
			div(Power);
			if (Pass == 2) 
				div(Power);
		} while( Zero != 0 );
	}
 
	if( (argc == 3) && (*argv[2] == 't') )
	{
		(void)fprintf(stderr,"end computation:   ");
		pitime();
	}
 
	if( printer )
	{
		(void)printf("\t\tThe value of pi" );
		(void)printf(" to %d decimal places\n\n\t\t",Size);
	}
 
	(void)printf("%d.",Result[0]);
	charcount = 0;
 
	for(i = 1; i <= Size; i++)
	{
		(void)printf( "%d", Result[i]);
		charcount++;
		if ( (charcount == 50) && printer )
		{
			(void)printf( "\n\t\t  " );
			charcount = 0;
		}
	}
	(void)printf( "\n" );
	return( 0 );	/* for lint */
}
 
add()
{
	register DIGIT *r, *t;
	register int sum, carry = 0;
    
	r = Result + Size;
	t = Term + Size;
 
	while( r >= Result )
	{
		sum = *r + *t + carry;
		carry = 0;
		if( sum >= 10 )
		{
			sum -= 10;
			carry++;
		}
		*r-- = sum;
		--t;
	}
}
 
sub()
{
	register DIGIT *r, *t;
	register int diff, loan = 0;		
	
	r = Result + Size;
	t = Term + Size;
 
	while( r >= Result )
	{
		diff = *r - *t - loan;
		loan =0;
		if( diff < 0 )
		{
			diff += 10;
			loan++;
		}
		*r-- = diff;
		--t;
	}
}
 
div(array)
register DIGIT *array;
{
	register DIGIT *end;
	register int quotient, residue, digit = 0;
 
	Zero = 0;
 
	for( end = array + Size; array <= end; )
	{
		digit    += *array;
		quotient =  digit / Divide;
		residue  =  digit % Divide;
 
		if((Zero !=0) || ((quotient+residue) != 0)) 
			Zero = 1;
		else 
			Zero = 0;
 
		*array++ = quotient;
		digit = 10 * residue;
	}
}
		
init()
{
	register DIGIT *p, *t, *r, *end;
 
	p = Power;
	t = Term;
	r = Result;
	end = Power+Size;
 
	while( p <= end )
	{
		*p++ = 0;
		*t++ = 0;
 
		if (Pass == 1) 
			*r++ = 0;
	}
 
	*Power = 16 / (Pass * Pass);
 
	if( Pass == 1 ) 
	{
		Divide = 5;
	}
	else
	{
		Divide = 239;
	}
 
	div(Power);
	Exp = 1;
	Sign = 3 - (2 * Pass);
}
 
copy()
{
	register DIGIT *t, *p, *end;
 
	t = Term;
	p = Power;
	end = Term + Size;
	
	while (t <= end)
	{
		*t++ = *p++;
	}
}
 
pitime()
{
	/* Print the current time \n */
 
	/* Here we just beep */
	(void)fprintf( stderr, "%c\n", 0x07 );
}