#define array_size 1000
#define TRUE 1
#define FALSE 0 

int IsPrime (int num);

                              /* shared data                      */
int a[array_size];            /*    array of numbers to sum       */
int global_index = 0;         /*    global index                  */
int long sum = 0;                  /*    global sum                    */
int long partial_sum = 0;          /*   local sum */

main ()
{
   long i;
   int nthreads,tid,work,mywork,leftover,start;
   int numProcs = 0;           /* Number of processors available  */
 
   omp_set_num_threads(3);

   numProcs = omp_get_num_procs();
   printf("There are %d processors currently available to the program.\n",numProcs);

   for (i = 0; i < array_size; i++)    // Initialize the array
      a[i] = i+1;

#pragma omp parallel default(shared) private (partial_sum,i,nthreads,tid,work,mywork,leftover,start)
{
  
   partial_sum = 0;
   tid = omp_get_thread_num();
   nthreads = omp_get_num_threads();
   work = array_size / (nthreads - 1);
   leftover = array_size % (nthreads - 1);

   if (tid == (nthreads - 1))
   {
      printf ("Number of threads = %d\n", nthreads);
      mywork = leftover;
   }
   else
      mywork = work;
   start = 0 + (tid * work);
   for (i = start; i < start + mywork ; i++)
   {  
    /* printf("Thread #%d adding a[%d] = %d to partial_sum of %d\n",
              tid, i, a[i],partial_sum);
    */
      if (IsPrime(a[i]))
         partial_sum += a[i];
   }

   #pragma omp critical
   {
      /* printf("Thread #%d adding partial_sum of %d to sum of %d\n",tid,partial_sum,sum); */
      sum += partial_sum;
   }
}

   printf("The sum of 1 to %i is %d\n", array_size, sum);
}  /* end of main */

int IsPrime (int num)
// precondition: num >=2
// postcondition: returns 1 if num is prime, else 0
{
    if (num == 2)         // the only even prime
        return TRUE;
    else if (num % 2 == 0)     // other even numbers are composite
        return FALSE;
    else
    {
        int prime = TRUE;
        int divisor = 3;
       // int upperLimit = static_cast(sqrt(num) + 1);
        int upperLimit = sqrt(num) + 1;
        while (divisor <= upperLimit)
        {
            if (num % divisor == 0)
                prime = FALSE;
            divisor +=2;
        }
        return prime;
    }
}