//
// The Pi Factory - RoB S 2013
//
// Computes Pi progressively using Machin's formula
//
// Columns in the output are:
//  DP = Decimal place within Pi
//  Pi  = The value of that decimal place
//  More Digits = Other stored digits which may not be accurate yet
//  Terms = Number of terms needed in each array so far
//
//  After calculation the final array contents are listed

#define ARRAYSIZE 1000    // Maximum number of terms (arbitrarily large here)
#define DECIMALPLACES 50  // Stops the programme 

#include <iostream.h>

int arccot1[ARRAYSIZE][3];    // In this version all arrays are the same size
//                        To save space they can be different sizes
//                        Actual sizes needed can be determined from the output
// In the array:
// Entry [n][0] contains the current remainder from y/(2n-1)
// Entry [n][1] contains the current remainder from 4*4*5/(5^2n)
// Entry [n][2] contains the current remainder from 4*239/(239^2n)

int piDigit[10]; // Last ten digits of pi so far - some may be inaccurate

main(){

	int n,m,l;  // loop counters
	int terms1=0,terms2=0,terms2a=0,dividend1,dividend2;
	int quotient1,quotient2,divisor2n,dividend2n,quotient2n;
	cout<<"The Pi Factory\n";
	cout<<"\t\t  Stored \t5^2n\t239^2n\n";
	cout<<"DP\tPi\t  Digits \tTerms\tTerms\n\n";

	// Initialise the arrays
	for(n=0;n<10;n++) piDigit[n]=0;
	for(n=0;n<ARRAYSIZE;n++){
		arccot1[n][0]=0;
		arccot1[n][1]=0;
		arccot1[n][2]=0;
		}

	// Calculate n times
	for (n=0;n<DECIMALPLACES+10;n++) {
		// Shift Pi along one digit
		for(m=0;m<9;m++) piDigit[m]=piDigit[m+1];
		piDigit[9]=0;
		// Scan all the terms and adjust Pi
		if(!n) dividend1=80; // First value is:
								 // 4 from pi/4 times 4 from Machin times 5 from Machin
		else dividend1=10*arccot1[1][1];
		if(!n) dividend2=956; // First value is:
								// 4 from pi/4 times 1 from Machin times 239 from Machin
		else dividend2=10*arccot1[1][2];
		// Process every term adding new ones if necessary
		for (m=1;m<=terms1||dividend1>0;m++ ){
			// Do the next step of division by 5^2
			quotient1=dividend1/25;
			// Calculate new /25 remainder for this term
			arccot1[m][1]=dividend1%25;
			// Calculate dividend for the next /25 division to use later
			dividend1=10*arccot1[m+1][1]+quotient1;

			if(m<=terms2||dividend2>0){
				// Do the next step of division by 239^2
				quotient2=dividend2/57121;
				// Calculate new remainder for this term
				arccot1[m][2]=dividend2%57121;
				// Calculate dividend for the next /57121 division to use later
				dividend2=10*arccot1[m+1][2]+quotient2;
				terms2a=m;
				quotient1-=quotient2;
				}

			// Calculate divisor for the 2m-1 division
			divisor2n=2*m-1;
			// Calculate dividend for the 2m-1 division
			// making it positive under all circumstances
			dividend2n=10*arccot1[m][0]+quotient1+9*divisor2n;

			// Calculate quotient for the 2m-1 division
			quotient2n=dividend2n/divisor2n-9;

			// Calculate new 2n-1 remainder for this term
			arccot1[m][0]=dividend2n%divisor2n;

			// Adjust pi for this term
				if(m%2==0)piDigit[9]-=quotient2n;   // Even terms subtract
				else piDigit[9]+=quotient2n;	       // and odd terms add

			// Do carries along the pi digit array - (a bit unsophisticated)
			for (l=9;l>0;l--){
				if(piDigit[l]<0){
					piDigit[l]+=10;
					piDigit[l-1]--;
					}
				else if(piDigit[l]>9){
					piDigit[l]-=10;
					piDigit[l-1]++;
					}
				else break;
				}
			// End of the terms in this array
			}
	  // New terms may have been added to the array, so update the terms counts
	  terms1=m-1;
	  terms2=terms2a;

	  // Print the current values
	  if(n>=9)cout<<(n-9);   // decimal place counter
	  else cout<<" ";
	  cout.width(1);
	  cout<<"\t"<<piDigit[0]
						<<"\t"<<piDigit[1]<<piDigit[2]<<piDigit[3]<<piDigit[4]
						<<piDigit[5]<<piDigit[6]<<piDigit[7]<<piDigit[8]
						<<piDigit[9];
	  cout.width(6);
	  cout<<"\t "<<terms1<<"\t "<<terms2<<"\n";

	  }
	  cout<<"\nSorry, no more pi's left.\n\n";
	  cout<<"Final remainders stored\n\n";
	  cout<<"Term\t5^2\t239^2\t2n-1\n";
	  cout<<"\tRem\tRem\tRem\n\n";
	  for(n=1;n<=terms1;n++){
		cout<<n<<"\t"<<arccot1[n][1]<<"\t"<<arccot1[n][2]<<"\t"<<arccot1[n][0]<<"\n";
		}

	  return 0;
}