/* 
   Compute integral from 0  to 1  4/(1+x*x). 
   This shoudl be equal to PI. 


   Parallel version using a parallel for construct. Note that we have
   critical race to the global sum, which we need to protect
   
   Scalability: good.

*/ 


#include <stdio.h> 
#include <stdlib.h> 
#include <omp.h>


#define NUM_THREADS 2

int main(int argc, char** argv) {

  long num_steps = 10000000; 
  double t1, t2; 
  
  double step = 1/(double)num_steps; 
  
  double global_sum = 0; 

  t1 = omp_get_wtime(); 
  
  int i;   
  omp_set_num_threads(NUM_THREADS); 
#pragma omp parallel  
  { 
    //private variables to each thread  
    double x, sum = 0; 

    int thread_id = omp_get_thread_num(); 
    printf("Start(%d)\n", thread_id); 
    
#pragma omp for 
    for (i=0; i<num_steps; i++) {
      //value in the middle of the rectangle
      x = (i+.5)*step; 
      sum = sum + 4/(1+x*x); 
    }
 
    //update global sum
#pragma omp critical 
    global_sum += sum; 

    printf("Done(%d)\n", thread_id); 
  }
  
  t2 = omp_get_wtime();   
  
  printf("Done. PI=%.3f, time = %.3f milliseconds\n", step*global_sum, (t2-t1)*1000);
  
  return 0; 
}