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


double f( double x ) 
{
  return 1/(2+sin(x)); 
}

double simpson( double a, double b )
{
  double h3, ab;
  h3 = (b-a)/6;
  ab = 0.5*(a+b);
  return h3*(f(a)+4*f(ab)+f(b));
}



main()
{
  int i, NUM_THREADS=3;
  double h, result;
  omp_set_num_threads(NUM_THREADS);
#pragma omp parallel reduction (+:result) 
  {
    //    NUM_THREADS = omp_get_num_threads();
    h = 3.0/NUM_THREADS;
    result=0.0;
    printf("I am number %d of %d\n",omp_get_thread_num(),NUM_THREADS);
#pragma omp for
    for ( i=0; i<NUM_THREADS; i++ ) 
      result += simpson(i*h,(i+1)*h);
  }
  printf("Approximate integral = %lf\n", result);
  return(0);
}