#include <stdio.h>
#include <math.h>
#define MAXLEVEL 30  /* Maximum level of recursion */

/* Adaptive and recursive Simpson's rule integration using
   Algorithm 145 by W. M. McKeeman, CACM 6, 167 (1963).
   Translated to C by H.L. Lynch 24 June 1983.

   Integrate the function "function" between lower_limit and upper_limit
   with a tolerance on the integration "tolerance".

   Calling sequence is 

      ans = integral (function, lower_limit, upper_limit,
                      tolerance, int_fail);

   Integral and all its arguments except int_fail are type double. The
   integer int_fail returns the number of times the recursion exceeded
   the maximum number allowed; 0 means that integration was successful.
   See tstintegral.c for an example of the use of integral. The function
   integral is just a nice user interface to simpson, which does all the
   work.  Simpson is a local function in the Algol original, but since C
   does not allow local functions, it appears as a static function here.
*/

/* Inter-function communication variables */
static int 
   level,       /* level of recursion reached */
   int_fails;   /* number of times recursion went too deep */

double integral (f, a, b, eps, int_fail)

   double (*f)();      /* function to be integrated */

   double a, b,        /* lower and upper limit of integration */ 
          eps;         /* tolerance on integration */

   int    *int_fail;   /* number of times recursion went too deep */

   /* Integral will numerically approximate the integral of the function
      f  between the limits a and b by the applicatiion of a modified
      Simpson's rule.  Although eps is a measure of the relative error
      of the result, the actual error may be very much larger (e.g. 
      whenever the answer is small because a positive area cancelled a
      negative area).  The procedure attempts to minimize the number of
      function evaluations by using small subdivisions of the interval
      only where required for the given tolerance.
   */

{  double simpson();
   double integral; /* Note: func name cannot be used without new dcl */

   level = 1;
   int_fails = 0;
   integral = simpson(f, a, b - a, (*f)(a),
                     4.0 * (*f)((a + b) / 2.0), (*f)(b), 1.0, 1.0, eps);
   /* Note weird mix of references to f in call to simpson. */
   *int_fail = int_fails;
   return (integral);
} /* end integral */

static double simpson(f, a, da, fa, fm, fb, absarea, est, eps)

   double (*f)();
   double a, da, fa, fm, fb, absarea, est, eps;

   /* This is the recursive Simpson's rule integration itself;
      simpson should be invisible outside this module.
   */

{  double dx, x1, x2, dx6, est1, est2, est3, f1, f2, f3, f4, sum;
   double integral, diff, tol, eps17, int1, int2, int3;

   dx = da / 3.0;
   x1 = a + dx;
   x2 = x1 + dx;
   f1 = 4.0 * (*f)(a + dx / 2.0);
   f2 = (*f)(x1);
   f3 = (*f)(x2);
   f4 = 4.0 * (*f)(a + 2.5 * dx);
   dx6 = dx / 6.0;
   est1 = (fa + f1 + f2) * dx6;
   est2 = (f2 + fm + f3) * dx6;
   est3 = (f3 + f4 + fb) * dx6;
   absarea = absarea - fabs(est) + fabs(est1) + fabs(est2) + fabs(est3);
   sum = est1 + est2 + est3;
   level = level + 1;
   /* Compiler bug: Following expression causes error during assembly: 
      if (((fabs(est - sum) <=  eps * absarea) &&
           (est !=  1.0)) || (level >=  20))
   */
   diff = fabs(est - sum);
   tol = eps * absarea;
   if (((diff <=  tol) && (est !=  1.0))  || (level >=  MAXLEVEL))
   {  integral = sum; 
      if (level >= MAXLEVEL) int_fails = int_fails + 1;
   }
   else
   {  /* Unless the integral is split up this way, the answer will go 
         crazy  if the recursion goes too deeply.  Reason not clear.
      */
      eps17 = eps / 1.7;
      int1     = simpson(f,  a, dx, fa, f1, f2, absarea, est1, eps17);
      int2     = simpson(f, x1, dx, f2, fm, f3, absarea, est2, eps17);
      int3     = simpson(f, x2, dx, f3, f4, fb, absarea, est3, eps17);
      integral = int1 + int2 + int3;
   }
   level = level - 1;
   return (integral); 
} /* end simpson */