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


   Parallel version using a parallel for construct and reduction. 
   
   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 sum = 0; 

  int i; 

  t1 = omp_get_wtime(); 
    
  omp_set_num_threads(NUM_THREADS); 
#pragma omp parallel  reduction(+:sum)
  { 
    //private variables to each thread  
    double x; 

#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); 
    }
  } //end parallel block 
  
  t2 = omp_get_wtime();   
  
  printf("Done. PI=%.3f, time = %.3f milliseconds\n", step*sum, (t2-t1)*1000);
  
  return 0; 
}