// Adaptive Quadrature with Simpsons Method

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

#define true 1
#define false 0

double func(double x)
{
   return (100.0 / ((x)*(x)) * sin(10.0/(x)));
}

void main()
{
   double tol[20], A[20], H[20], FA[20], FC[20], FB[20], S[20], V[7];
   int L[20];
   double APP=0.0, FD, FE, S1, S2;
   int count=0, I=1, LEV, OK=true;
   double left_end=1.0, right_end=1.5; // 좌우측 경계치
   double eps=1.0e-6;     // 오차한계
   double num=20;         // maximum number of levels

   if (OK) {
      tol[I] = 10.0 * eps;
      A[I] = left_end;
      H[I] = 0.5 * (right_end - left_end);
		FA[I] = func(left_end), count++;
      FC[I] = func(left_end+H[I]), count++;
		FB[I] = func(right_end), count++;

      /* Approximation from Simpsons method for entire interval. */
      S[I] = H[I] * (FA[I] + 4.0 * FC[I] + FB[I]) / 3.0;
      L[I] = 1;

      while ((I > 0) && OK) {
    	   FD = func((A[I] + 0.5 * H[I])), count++;
	   FE = func((A[I] + 1.5 * H[I])), count++;
	   /* Approximation from Simpsons method for halves of intervals.*/
	   S1 = H[I] * (FA[I] + 4.0 * FD + FC[I]) / 6.0;
	   S2 = H[I] * (FC[I] + 4.0 * FE + FB[I]) / 6.0;
	   /* Save data at this level */
	   V[0] = A[I];
	   V[1] = FA[I];
	   V[2] = FC[I];
	   V[3] = FB[I];
	   V[4] = H[I];
	   V[5] = tol[I];
	   V[6] = S[I];
	   LEV = L[I];

	   /* Delete the level */
	   I--;

	   if (fabs(S1 + S2 - V[6]) < V[5]) APP = APP + (S1 + S2);
	   else {
	       if (LEV >= num) OK = false;   /* procedure fails */
	    	  else {
	           /* Add one level */
	           /* Data for right half subinterval */
	           I++;
	           A[I] = V[0] + V[4];
	           FA[I] = V[2];
	           FC[I] = FE;
	           FB[I] = V[3];
	           H[I] = 0.5 * V[4];
	           tol[I] = 0.5 * V[5];
	           S[I] = S2;
	           L[I] = LEV + 1;

	           /* Data for left half subinterval */
	           I++;
	           A[I] = V[0];
	           FA[I] = V[1];
	           FC[I] = FD;
	           FB[I] = V[2];
	           H[I] = H[I-1];
	           tol[I] = tol[I-1];
	           S[I] = S1;
		   L[I] = L[I-1];
	      }
	   }
      }
      if (!OK) printf("Level exceeded\n");
	else {
	   printf("\nThe integral of F from %12.8f to %12.8f is\n", left_end, right_end);
   	   printf("%12.8f to within %14.8e\n", APP, eps);
	   printf("The number of function evaluations is: %d\n", count);
      }
   }
}