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

#define MAX 32000


double simpson(double y[],double, double, int);

int main()
{

  int i,M;
  double a,b,s;
  double x[MAX],y[MAX];

  M=10;
  a=1.0;
  b=6.0;

  for(i=0;i<=2*M;i++)
    {
      x[i]=a+(b-a)*i/(2.0*M);
      y[i]=2.0+sin(2.0*sqrt(x[i]));
    }

  s=simpson(y,a,b,M);

  printf("%f\n",s);

  return 1;

}


double simpson(double y[],double a,double b,int M)
{

  /* 
     Calculo de la integral por la regla compuesta
     de Simpson para la funcion y[0]...y[2*M] y
     a=x[0] b=x[2*M].

     El error es del orden de [(b-a)/(2*M)]^4
 
  */

  int k;
  double h,s1,s2,s;

  s1=0.0;
  s2=0.0;

  h=(b-a)/(2.0*M);

  for(k=1;k<=M;k++)
    {
      s1=s1+y[2*k-1];
    }

  for(k=1;k<=M-1;k++)
    {
      s2=s2+y[2*k];
    }

  s=h*(y[0]+y[2*M]+4*s1+2*s2)/3.0;

  return s;

}