/*
   Name: Plouffe and Bellard algorithm, version 1
   Author: Fabrice Bellard, http://fabrice.bellard.free.fr/pi/
   Cosmetic Revision: Sayan Chakraborti, http://sourceforge.net/users/sayanchak/
   Description: This program calculates decimal digits of Pi
   Date: Monday, May 13, 2003
   Copyright: This program is distributed under the GNU General Public License

   OpenMosix addition
   Author: Jeremy Perrin, picalcman@yahoo.com
   Description: Adds OpenMosix capablities to above PI calculator
   Date: Thursday, September 28, 2003
   Copyright: This program is distributed under the GNU General Public License
   Coming Soon: Threaded OpenMosix capable version, using FSU_Thread library
*/

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

#include <sys/wait.h> //Needed for waitpid

//Includes for Cluster Node Counting
#include <sys/types.h>
#include <string.h>
#include <unistd.h>
#include <dirent.h>
#define clusterdir "/proc/hpc/nodes"

#ifdef HAS_LONG_LONG
#define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
#else
#define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
#endif

typedef struct{
	int PID;	//Process Identification
	int StatFlag;	//Return status of PID from waitpid
} PIDSaver;

/*///////////////////////////////////////////////////////////////
Returns the number of CPU's in the OpenMosix Cluster
Fucntion source provided by
http://www.openmosixview.com/docs/openMosixAPI.html
*////////////////////////////////////////////////////////////////
int CPUCount() {
	int howmany=0;
	int cpus=0;

	char tmpdirname[200];
	char tmpfilename[200];

	struct dirent *dir_info;

	FILE *fp;
	DIR *dirhpc;

	strcpy(tmpdirname, "");

	if ((dirhpc=opendir(clusterdir))!=NULL) {
		while ((dir_info = readdir(dirhpc))!=NULL) {

			strcpy(tmpdirname, dir_info->d_name);
			strcpy(tmpfilename, "/proc/hpc/nodes/");

			strcat(tmpfilename, tmpdirname);
			strcat(tmpfilename, "/cpus");

			if (!strchr(tmpdirname, '.')) {
				fp=fopen(tmpfilename, "r");

				if (fp) {
					fscanf(fp, "%d", &cpus);

					if(cpus>0) {
						howmany=howmany+cpus;
					}
					fclose(fp);
				}
			}
		}
		closedir(dirhpc);
		printf("OpenMosix directory found, Node count is %d\n", howmany);

		return howmany;
	}
	else {
		printf("Can't find OpenMosix directory: %s\nForcing Node count to 2\n", clusterdir);
		return(2);
	}
}


/* returns the inverse of x mod y */
int inv_mod(int x,int y){
	int q,u,v,a,c,t;

	u=x;
	v=y;
	c=1;
	a=0;

	do {
		q=v/u;

		t=c;
		c=a-q*c;
		a=t;

		t=u;
		u=v-q*u;
		v=t;
	} while (u!=0);

	a=a%y;
	if (a<0) a=y+a;
	return a;
}

/* returns (a^b) mod m */
int pow_mod(int a,int b,int m){
	int r,aa;

	r=1;
	aa=a;

	while (1) {
		if (b&1) r= (int)mul_mod(r,aa,m);
		b=b>>1;
		if (b == 0) break;
		aa= (int)mul_mod(aa,aa,m);
	}

	return r;
}

/* returns true if n is prime */
int is_prime(int n){
	int r,i;

	if ((n % 2) == 0) return 0;

	r=(int)(sqrt(n));
	for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;

	return 1;
}

/* returns the prime number immediatly after n */
int next_prime(int n)
{
	do {
		n++;
	} while (!is_prime(n));

	return n;
}

/* finds out nine digits of pi starting from n */
/* will probably be the main function in openmosix capable version*/
int pseudo_main(int n){

	int av,a,vmax,N,num,den,k,kq,kq2,t,v,s,i;
	double sum;
	char filename[100];
	FILE *f_out;

	double start_time, elap_time;
	start_time = (double) clock ();

	sprintf (filename, "%i.tmp", n);
	f_out = fopen (filename, "w");

	N=(int)((n+20)*log(10)/log(2));

	sum=0;

	for(a=3;a<=(2*N);a=next_prime(a)) {

		vmax=(int)(log(2*N)/log(a));
		av=1;

		for(i=0;i<vmax;i++) av=av*a;

		s=0;
		num=1;
		den=1;
		v=0;
		kq=1;
		kq2=1;

		for(k=1;k<=N;k++) {

			t=k;
			if (kq >= a) {
				do {
					t=t/a;
					v--;
				} while ((t % a) == 0);
				kq=0;
			}

			kq++;
			num = (int)mul_mod(num,t,av);

			t=(2*k-1);
			if (kq2 >= a) {
				if (kq2 == a) {
					do {
						t=t/a;
						v++;
					} while ((t % a) == 0);
				}
				kq2-=a;
			}

			den = (int)mul_mod(den,t,av);
			kq2+=2;

			if (v > 0) {
				t=inv_mod(den,av);
				t= (int)mul_mod(t,num,av);

				t= (int)mul_mod(t,k,av);
				for(i=v;i<vmax;i++) t= (int)mul_mod(t,a,av);
				s+=t;
				if (s>=av) s-=av;
			}

		}

		t=pow_mod(10,n-1,av);
		s= (int)mul_mod(s,t,av);
		sum=fmod(sum+(double) s/ (double) av,1.0);
	}

	elap_time = ((double) clock () - (double) start_time) / CLOCKS_PER_SEC;

	printf("Position: %9d Digits: %09d Time: %9.2lf seconds\n",n,(int)(sum*1e9),elap_time);
	fprintf (f_out, "%09d\n", (int)(sum*1e9));
	fclose (f_out);

	return(0);
}



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

	int position=1,counter=1,number;
	int NodeCount =0, ProcessSpawned =0;  //Added by Jeremy Perrin, track Number of nodes
	double start_time, elap_time;
	char filename_in[100], filename_out[100];
	FILE *f_out, *f_in;

	if (argc<2 || (number=atoi(argv[1])) <= 8){
		printf("\nPlease type the number of digits of pi you want, after the filename.\n");
		printf("Since you have not done that please type it now.\nNumber of Digits = ");
		scanf("%d",&number);
	}

	start_time = (double) clock();

	number = number / 9;
	number = number * 9;

	printf("\nComputing %d digits of pi ...\n\n",number);

	/*
	Child Spawning process algorithm by Jeremy Perrin, picalcman@yahoo.com
	*/

	//Find number of OpenMosix Nodes
	NodeCount = CPUCount();
	NodeCount *= 6; //Six processes per Node

	//Make sure there isn't more processes than pi calculations 
	if(NodeCount > (number/9)){
		NodeCount = number / 9;
		printf("Too many processes, cutting back to %d processes\n", NodeCount);
	}

	//Create dynamic sized array
	PIDSaver* childpids = (PIDSaver*)new PIDSaver[NodeCount];

	//Initialize dynamic array componets
	printf("Initializing Process Structures: ");
	for(int j=0; j < NodeCount; j++){
		childpids[j].PID = 0;
		childpids[j].StatFlag = 0;
		printf("+");
	}
	printf("\n\n");

	counter = 1;
	while(counter < number){
		for(int j =0; j < NodeCount; j++){
			//Check if child process is running if not start a new one, save PID
			int p = waitpid(childpids[j].PID, &childpids[j].StatFlag, WNOHANG);
			if(p == childpids[j].PID || p < 0)
			{
				if(counter < number){
					//Fork new process with pseudo_main, save PID
					if((childpids[j].PID = fork()) == 0){
						pseudo_main(counter);
						exit(0);
					}
					ProcessSpawned++; //Count spawned processes
					counter += 9;
				}
				else{
					j = NodeCount; //Force exit of for loop
				}
			}
		}
	}

	//Wait for all child process to end
	for(int k =0; k < NodeCount; k++){
		waitpid(childpids[k].PID, &childpids[k].StatFlag, 0);
	}

	//Done with childpids, delete dynamic array
	delete childpids;

	printf("\nTotal spawned processes: %d\n", ProcessSpawned);

	/*
	End of child Spawning algorithm by Jeremy Perrin, picalcman@yahoo.com
	*/

	elap_time = ((double) clock () - (double) start_time) / CLOCKS_PER_SEC;

	printf("\nTotal calculation done in: %9.2lf seconds\n",elap_time);

	printf("\nGenerating output file ...\n");
	start_time = (double) clock();

	sprintf (filename_out, "pi%i.txt", number);
	f_out = fopen (filename_out, "w"); fprintf(f_out, "3.");

	counter = 1;
	while (counter<number){
		sprintf (filename_in, "%i.tmp", counter);
		f_in = fopen (filename_in, "r");
		
		position = 1;
		while (position<10){
			putc((getc (f_in)),f_out);
			position=position+1;
		}
		
		counter=counter+9;
		fclose (f_in);

/* Next line contributed by Arko Provo Mukherjee, arkoprovo@peacefulaction.org */
		remove (filename_in);
	}

	fclose (f_out);

	elap_time = ((double) clock () - (double) start_time) / CLOCKS_PER_SEC;
	printf("\n%s generated in: %9.2lf seconds\n\n",filename_out,elap_time);
	
	exit(0);
}