//Cuadratura adaptativa basada en la Regla de Simpsom
//para el calculo de integrales
//a, b, son respectivamente el extremo inferior
//y el superior del intervalo de integracion
//N es el numero de intervalos
//TOL, es la presicion
//Ramiro alcocer
//valcoey@hotmail.com
//www.oocities.org/valcoey/index.html

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

//maximo numero de pasos  (intervalos)
#define N	100
//valor de error
#define TOL	0.0000001

double f(double x)
{
 return (exp(-x*x));
}

int main(void)
{
int i;
double a, b, sum;
double s1, s2;
double v1, v2, v3, v4, v5, v6, v7, v8;
double fe, fd;
double aa[N], h[N], fa[N], fc[N], fb[N];
double s[N], L[N], tol[N];

clrscr();
//Extremo inferior del intervalo
a = 1.0;
//Extremo superior del intervalo
b = 1.5;
sum=0.0;
i=1;
tol[i]=10.0*TOL;
aa[i]=a;
h[i]=(b-a)/2.0;
fa[i]=f(a);
fc[i]=f(a+h[i]);
fb[i]=f(b);
s[i]=h[i]*(fa[i] + 4.0*fc[i] + fb[i])/3.0;
L[i]=1.0;
while (i>0)
	{
	fd=f(aa[i] + h[i]/2.0);
	fe=f(aa[i] + 3.0*h[i]/2.0);
	s1=h[i]*(fa[i] + 4.0*fd + fc[i])/6.0;
	s2=h[i]*(fc[i] + 4.0*fe + fb[i])/6.0;
	v1=aa[i];
	v2=fa[i];
	v3=fc[i];
	v4=fb[i];
	v5=h[i];
	v6=tol[i];
	v7=s[i];
	v8=L[i];
	i=i-1;
	if (fabs(s1+s2-v7)<v6)
		sum=sum + (s1 + s2);
	else
		if (v8>=N)
			{
			printf("ERROR: Numero insuficiente de intervalos");
			getch();
			return 0;
			}
		else
			{
			i=i+1;
			aa[i]=v1+v5;
			fa[i]=v3;
			fc[i]=fe;
			fb[i]=v4;
			h[i]=v5/2.0;
			tol[i]=v6/2.0;
			s[i]=s2;
			L[i]=v8+1;
			i=i+1;
			aa[i]=v1;
			fa[i]=v2;
			fc[i]=fd;
			fb[i]=v3;
			h[i]=h[i-1];
			tol[i]=tol[i-1];
			s[i]=s1;
			L[i]=L[i-1];
			}
		}
//valor de la integral
printf("%.17lf\n", sum);
getch();
return 0;
}