/*
 *	Sieve of Eratosthenes
 *
 *	Programmed by Michael J. Quinn
 *
 *	Last modification: 7 September 2001
 *
 *	Modified by Tomofumi Yuki
 *
 */

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MILLION 1000000L

void usuage(char*);

int main (int argc, char *argv[]) {
	int	 count;				/* Local prime count */
	double elapsed_time;	/* Parallel execution time */
	int  first;				/* Index of first multiple */
	int  global_count;		/* Global prime count */
	int  high_value;		/* Highest value on this proc */
	int  i;
	int  id;				/* Process ID number */
	int  index;				/* Index of current prime */
	int  low_value;			/* Lowest value on this proc */
	char *marked;			/* Portion of 2,...,'n' */
	int  n;					/* Sieving from 2, ..., 'n' */
	int  p;					/* Number of processes */
	int  proc0_size;		/* Size of proc 0's subarray */
	int  prime;				/* Current prime */
	int  size;				/* Elements in 'marked' */

	int *fprimes;			/* For Output first 3 primes*/
	int *lprimes;			/* For Output last 3 primes*/
	int j;

	struct timeval start;
	struct timeval end;

	/* Initialize MPI */
	if (MPI_Init (&argc, &argv) || 
		MPI_Comm_rank (MPI_COMM_WORLD, &id) || 
		MPI_Comm_size (MPI_COMM_WORLD, &p)) {
		if (!id) printf("MPI Initialize failed.\n");
		MPI_Finalize();
		exit(1);
	}

	/* Start the timer */
	MPI_Barrier(MPI_COMM_WORLD);
	gettimeofday(&start, NULL);

	/* Input Check */
	if (argc != 2) {
		if (!id) usuage(argv[0]);
		MPI_Finalize();
		exit (1);
	}

	n = atoi(argv[1]);
	if (n < 2) {
		printf ("n must be greater than 1.\n");
		exit(1);
	}

	/* Figure out this process's share of the array, as
		well as the integers represented by the first and
		last array elements */
	size = n/p;

	low_value = size*id;
	high_value = size*(id+1);

	/*Last processor needs to take the remainder */
	if (id==p-1) {
		high_value=n;
		size = high_value - low_value + 1;
	}


	/* Bail out if all the primes used for sieving are
		not all held by process 0 */
	proc0_size = (n-1)/p;

	if ((2 + proc0_size) < (int) sqrt((double) n)) {
		if (!id) printf ("Too many processes.\n");
		MPI_Finalize();
		exit (1);
	}

	/* Allocate this process's share of the array. */

	marked = (char *) malloc (size);

	if (marked == NULL) {
		printf ("Cannot allocate enough memory.\n");
		MPI_Finalize();
		exit (1);
	}

	/*Initialize array marked*/
	for (i = 0; i < size; i++) marked[i] = 0;
	if (!id) index = 0;
	if (!id) marked[0] = marked[1] = 1; /*0 and 1 starts marked*/
	prime = 2;

	do {
		if (prime * prime > low_value)
			first = prime * prime - low_value;
		else {
			if (!(low_value % prime)) first = 0;
			else first = prime - (low_value % prime);
		}
		/*Marks off prime*x*/
		for (i = first; i < size; i += prime) marked[i] = 1;
		if (!id) {
			while (marked[++index]);
			/*Next primes is the non marked index*/
			prime = index;
		}
		if (p > 1) MPI_Bcast (&prime,  1, MPI_INT, 0, MPI_COMM_WORLD);
	} while (prime * prime <= n);

	/*Count the primes for each process*/
	count = 0;
	for (i = 0; i < size; i++) {
		if (!marked[i]) count++;
	}

	/*Collect data from all processess*/
	if (p > 1) {
		MPI_Reduce (&count, &global_count, 1, MPI_INT, MPI_SUM,
	  		 0, MPI_COMM_WORLD);
	} else {
			  global_count = count;
	}

	/* Stop the timer */
	gettimeofday(&end, NULL);

	elapsed_time = (end.tv_sec - start.tv_sec + (double)(end.tv_usec - start.tv_usec)/MILLION);

	/* Print the results */

	#ifdef DEBUG
	/*Get first and last three primes*/
	fprimes = (int *)malloc(sizeof(int) * 3);
	lprimes = (int *)malloc(sizeof(int) * 3);

	for (i = 0; i < 3; i++) {
		fprimes[i] = 0;
		lprimes[i] = 0;
	}

	j = 0;
	for (i = 0; i < size; i++) {
		if (!marked[i]) {
			fprimes[j] = i + low_value;
			j++;
		}
		if (j > 2) {
			break;
		}
	}

	j = 0;
	for (i = size - 1; i >= 0; i--) {
		if (!marked[i]) {
			lprimes[j] = i + low_value;
			j++;
		}
		if (j > 2) {
			break;
		}

	}

	printf ("%d: %d %d %d %d %d %d\n", id, 
		fprimes[0], fprimes[1], fprimes[2],
		lprimes[2], lprimes[1], lprimes[0]);
	MPI_Barrier(MPI_COMM_WORLD);

	if (!id) {
		printf ("There are %d primes less than or equal to %d\n",
			global_count, n);
	}
	free(fprimes);
	free(lprimes);
	#endif

	if (!id) {
       		printf("Elapsed Time=%lf\n", elapsed_time);
	}

	free(marked);
	MPI_Finalize ();
	return 0;
}

void usuage(char *p_name) {
	printf ("Usuage: %s n\n\n", p_name);
	printf ("Calculates the number of prime numbers up to n.\n");
}