#include <stdio.h>
#include <mpi.h>
#define N 100
#define TRUE 1
#define FALSE 0
#define nonprime 0
#define Square(x) ((x)*(x))

void main(int argc, char **argv )
 {
  int *a, *aprime, *primes;
  int i, k, p, my_k;
  int found, myid, numprocs, myfirst, mylast;
  double t1, t2;
  MPI_Status status;

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid); 

  // 1. create a list of natural numbers 2, 3, 4, ... none of which is marked.
  a = (int *) malloc (N * sizeof(int));
  aprime = (int *) malloc (N * sizeof(int));
  myfirst = myid*N;
  for (i=0;i<N;i++) { a[i] = myfirst+i; aprime[i] = a[i]; } 
  // master marks zero & one as non-prime
  if (myid == 0)
   {
    aprime[0] = nonprime;
    aprime[1] = nonprime;
    my_k = a[2]; 
   }
  else my_k = a[0];
  k = 1;

  t1 = MPI_Wtime();
  // repeat while k^2 < N
  while (Square(k) <= (N*numprocs))
   {
    // find the mimimum k from all processes
    MPI_Allreduce(&my_k,&k,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
    // a. Mark all multiples of k between k^2 and N
    myfirst = 0;
    mylast = N;
    if ((Square(k) >= a[0]) && (Square(k) <= a[N-1]))// Square(k) is in my range
     myfirst = Square(k) % N;
    if (Square(k) > a[N-1])// Square(k) is greater than my range
      myfirst = mylast;

    for (i=myfirst; i<mylast; i++) // serial: for (i=Square(k);i<N;i++)
     if (aprime[i] && ((a[i] % k) == 0)) { aprime[i] = nonprime; }
    // b. Find the smallest number greater than k that is unmarked 
    //    and set k to this new value until k^2 > N
    found = FALSE;
    for (i=0;!found;i++)
     {
      if (aprime[i] && (a[i] > k))
       {
        my_k = a[i]; 
        found = TRUE; 
       }
     }
   }
  t2 = MPI_Wtime();
  // do a gather & have the master print out all primes
  if (myid == 0) primes = (int *) malloc(numprocs * N * sizeof(int));
  MPI_Gather(aprime,N,MPI_INT,primes,N,MPI_INT,0,MPI_COMM_WORLD);
  // 4. The unmarked numbers are primes
  p = 0;
  if (myid == 0)
   {
    for (i=0;i<(N*numprocs);i++) 
     if (primes[i]) { printf("%d\t%d\n",i,primes[i]); p++; } 
    printf("%d primes found in %.2f seconds\n",p,t2-t1);  
   }
  MPI_Finalize();
 }