#include <pthread.h>
#include <stdio.h>
#include <math.h>
#define NUM_THREADS 3

pthread_mutex_t reduction_mutex;
pthread_t *tid;
double result=0.0;

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));
}

void *worker( void *arg )
{
  int my_id;
  double my_result;

  my_id = *((int*)arg);
  my_result = simpson( (double)my_id, (double)my_id+1.0 );
  printf("Worker #%d got %lf \n", my_id, my_result);

  pthread_mutex_lock( &reduction_mutex );
  result += my_result;
  pthread_mutex_unlock( &reduction_mutex );
}

main()
{
  int i;
  pthread_t tid[NUM_THREADS];
  int t_num[NUM_THREADS]={0,1,2};

  pthread_mutex_init( &reduction_mutex, NULL );

  for ( i=0; i<NUM_THREADS; i++ ) 
    pthread_create( &tid[i], NULL, worker, (void*)&t_num[i] );

  for ( i=0; i<NUM_THREADS; i++ )
    pthread_join( tid[i], NULL );

  printf("Approximate integral = %lf\n", result);
  return(0);
}