#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define TYPE double
#define FUNC func

TYPE midpoint(TYPE a, TYPE b, int n);
TYPE midpoint2(TYPE a, TYPE b, int n); 
double mcIntegration(double a, double b, int N);
TYPE func(TYPE);

int main()
{
  srand48( time(0) ); 

  TYPE a=0., b=1.;
  int M= 10;
  int i;

  /* for(i=1;i<20;i++) */
  /*   { */
  /*     M=i; */
  /*     TYPE I=midpoint(a,b,M); */
  /*     // printf("I: %lf \n",I); */
  /*     TYPE I2=midpoint2(a,b,M); */
  /*     // printf("I2: %lf \n",I2); */
  /*     printf("%d \t %lf \t %lf \n",M,I,I2); */
  /*   } */

  /* for(i=1;i<1000;i++) */
  /*   { */
  /*     M=i*100; */
  /*     TYPE I3=mcIntegration(a,b,M); */
  /*     printf("%d \t %lf \n", M, I3); */
  /*   } */

  TYPE I=midpoint(a,b,M);
  printf("I: %lf \n",I);
  TYPE I2=midpoint2(a,b,M);
  printf("I2: %lf \n",I2);

  /* TYPE I3=mcIntegration(a,b,M*10); */
  /* printf("I3: %lf \n",I3); */
  

  /* float a1=0.1; */
  /* double a2=0.1; */
  /* printf("diff: %e \n", a1-a2); */

  return 0;
}

TYPE midpoint(TYPE a, TYPE b, int n) 
{
  /* una cattiva realizzazione dell'algoritmo del punto di mezzo */
  TYPE I = 0., x = a;
  TYPE h = (b - a) / (TYPE) n;

  int stepcounter=0;

  while (x < b) 
    { 
      I += FUNC(x + 0.5 * h);
      //
      x += h;
      stepcounter++;
    }

#ifdef _DEBUG
  printf("midpoint h: %.55f, n steps: %d  \n", h, stepcounter );  
#endif
  I *= h;
  return  I;
}


TYPE midpoint2(TYPE a, TYPE b, int n) 
{
  /* Una buona realizzazione dell'algoritmo del punto di mezzo */
  TYPE I = 0., x; 
  TYPE h = (b - a) / (TYPE) n;
  int i;

  int stepcounter=0;

  for (i = 0; i < n; i++) 
    {
      x = a + h * i + 0.5 * h;
      I += FUNC(x);
      stepcounter++;
    }

#ifdef _DEBUG
  printf("midpoint2 h: %.55f n steps: %d \n", h, stepcounter );
#endif

  return I * h;
}
 
double mcIntegration(double a, double b, int N) 
{
  int i;
  double S = 0.;
  double x;
  
  for (i = 0; i < N; i++) 
    {
      x = a + (b - a) * (double)lrand48() / (double)RAND_MAX;
      S += FUNC(x);
    }
  return S * (b - a) / N;
}


TYPE func(TYPE x)
{
  return sin(x);
  //return x*x;
}