#include <stdio.h>
#include <omp.h>
#include <math.h>
#define N 100000000
#define TRUE 1
#define FALSE 0

int main(int argc, char **argv )
{
  char host[80];
  int *a, *p;
  int i, j, k, index_max, threads, pcount;
  double t1, t2;
  int found;
  int index_sieve;
  int jlimit, index_clump=2200000, ii, iistart, iistop;

  /* Set number of threads equal to argv[1] if present */
  if (argc > 1)
    {
      threads = atoi(argv[1]);
      if (omp_get_dynamic())
        {
          omp_set_dynamic(0);
          printf("called omp_set_dynamic(0)\n");
        }
      omp_set_num_threads(threads);
    }
  printf("%d threads max\n",omp_get_max_threads());

  index_sieve = (sqrt(N)+1)/2;
  index_max = (N-1)/2;

  a = (int *) malloc((index_max + 1) * sizeof(int));
  // 1. Create a list of odd numbers 3, 5, 7, ... none of which is marked.
  //    Index j corresponds to the odd number q = 2*j + 1.
  //    Conversely, odd number q corresponds to index j = (q-1)/2.
  //    If N is even, (N-1)/2 truncates to (N-2)/2 = highest odd number < N.
  for (j=1;j<=index_max;j++) a[j] = 1;

  // 2. Set k = 3, the first unmarked number on the list.
  k = 3;

  t1 = omp_get_wtime();
  // 3. Repeat
  //    FIRST PASS: find all primes less than sqrt(N)
  #pragma omp parallel firstprivate(k) private(i,found)
  while (k*k <= N)
    {
      // a. Mark all ODD multiples of k between k^2 and N;
      //    stride 2*k in natural numbers corresponds to stride k in index
      #pragma omp for
      for (i=(k*k-1)/2; i<index_sieve; i+=k) a[i] = 0;
      // 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=(k+1)/2;!found;i+=1)
      {
        if (a[i]){ k = 2*i+1; found = TRUE; }
      }
    }
  t2 = omp_get_wtime(); 
  printf("%.6f seconds\n",t2-t1);

  // 4. The unmarked numbers are primes
  pcount = 1;
  p = (int *) malloc(sqrt(N)*sizeof(int));
  //  printf("%d\n",2);
  for (i=1;i<index_sieve;i++)
    {
      if( a[i] )
        { 
          pcount++;
          p[pcount-2] = 2*i+1;
          //printf("%d\n",2*i+1);
        }
    }

  jlimit=pcount-1;
  t1 = omp_get_wtime();
  // 3. Repeat
  //    SECOND PASS: sieve cache-size clumps, given all the primes < sqrt(N)
  for (i=index_sieve; i<=index_max; i+=index_clump)
    {
      for (j=0;j<jlimit;j++)
        { 
          iistart = (2*i+1)/p[j];
          iistart = ((iistart + 1 - iistart%2)*p[j] - 1)/2;
          iistop = i + index_clump;
          if (iistop >= index_max) iistop = index_max + 1;
          #pragma omp parallel for
          for (ii=iistart; ii<=iistop; ii+=p[j]) a[ii]=0;
        }
    }
  t2 = omp_get_wtime(); 
  printf("%.6f seconds\n",t2-t1);

  // 4. The unmarked numbers in the rest of the list are also primes
  for (i=index_sieve; i<=index_max; i++)
    {
      if( a[i] )
        { 
          pcount++;
          //printf("%d\n",2*i+1);
        }
    }
  printf("%d primes between 0 and %d\n",pcount,N);
}