/*
The algorithm is very simple:
For a given n, the numbers <= n is not a prime if it can be factor by a 
number, say p, and p must be less than or equal to sqrt of n.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>

#define n 1000000

int main(int argc, char *argv[])
{
	long i, j, xs, ps, ms, Lbound = 0, Ubound = 0, temp;
	double clock0, clock1;
	float time;
	long m, s, count, total;
	long *p, *x;

	int rank, processes;

	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &processes); /* Getting # of processes */
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);     /* Getting my rank        */

	MPI_Barrier(MPI_COMM_WORLD); /* Make sure processes are synced */
	clock0=MPI_Wtime();

	m = floor(sqrt(n));
	s = floor((m-3)/2) + 1;

	ps = s;			/*size of p*/

	/* initialization */
	p = calloc(s, sizeof(long));
	for (i = 0; i < s; i ++)
		p[i] = 0;

	j = 0;
	for (i = 3; i <= m; i += 2)
	{
		p[j] = i;
		j ++;
	}

	ms = p[ps - 1];	/*ms = largest element in p*/

	s = floor((n-ms-2)/2)+1;

	xs = s;			/*size of x*/

	/* initialization */
	x = calloc(s, sizeof(long));
	for (i = 0; i < s; i ++)
		x[i] = 0;

	j = 0;
	for (i = ms + 2; i <= n; i += 2)
	{
		x[j] = i;
		j ++;
	}

	/* calculation */
	for (i = 0; i < ps - 1; i++)
	{
		for (j = i + 1; j < ps; j++)
		{
			if (p[i] != 0 && p[j] % p[i] == 0)
			{
				p[j] = 0;
			}
		}
	}

	temp = ceil(xs / processes);
	if (rank == 0)
		Lbound = 0;
	else
		Lbound = rank * temp;

	if (rank == processes - 1)
		Ubound = xs;
	else
		Ubound = (rank + 1) * temp;
/*
	if (rank == 0)
	{
		printf("xs %d\n", xs);
	}
	printf("%d %d %d\n", rank, Lbound, Ubound);
	(void) fflush(stdout);
*/
	for (i = 0; i < ps; i++)
	{
		for (j = Lbound; j < Ubound; j++)
		{
			if (p[i] != 0 && x[j] % p[i] == 0)
			{
				x[j] = 0;
			}
		}
	}
	MPI_Barrier(MPI_COMM_WORLD); /* Make sure processes are synced */
	clock1=MPI_Wtime();

	/* prepare the results for display */
	count = 0;
	if (rank == 0)
	{
		for (i = 0; i < ps; i ++)
		{
			if (p[i] != 0)
			{
/*				printf("p[%d]:%d\n", i, p[i]);
*/				count++;
			}
		}
	}

	for (i = Lbound; i < Ubound; i ++)
	{
		if (x[i] != 0)
		{
/*			printf("x[%d]:%d\n", i, x[i]);
*/			count++;
		}
	}

	MPI_Reduce(&count, &total, 1, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD);

	if (rank == 0)
	{
		/* display the result */
		printf("n :%d\n", n);
		printf("no. of prime :%d\n", total + 1);
		printf("time used: %e sec.\n", clock1 - clock0);
	}
	MPI_Finalize();

	free(p);
	free(x);
	return(0);
}