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

int main(int argc, char **argv )
{
  char host[80];
  int *a;
  int i, j, k, index_max, threads, pcount;
  double t1, t2;
  int found;

  /* 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_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
  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 parallel for
      for (i=(k*k-1)/2; i<=index_max; 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;
  //  printf("%d\n",2);
  for (i=1;i<=index_max;i++)
    {
      if( a[i] )
        { 
          pcount++;
          //printf("%d\n",2*i+1);
        }
    }
  printf("%d primes between 0 and %d\n",pcount,N);
}