function integral (function f(x : real) : real;  a, b, eps : real;
                   var int_fail : integer
                  ) : real;

   {  Adaptive and recursive Simpson's rule integration using
      Algorithm 145 by W. M. McKeeman, CACM 6, 167 (1963).
      Translated to Pascal by H.L. Lynch 24 March 1984.

      Integral will numerically approximate the integral of a
      function_f between the limits upper_limit and lower_limit
      by the applicatiion of a modified Simpson's rule.
      Although tolerance 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.

      Calling sequence is

         ans := integral(function_f, lower_limit, upper_limit,
                         tolerance, failures);

      failures is returned 0 if the maximum level of recursion was
      never reached in the integration; otherwise it returns the number
      of times the maximum recursion level was reached attempting to
      satisfy the tolerance.

      The function integral is just a nice user interface to simpson,
      which does all the work.
   }

const
   MAXLEVEL = 30;  { Maximum level of recursion }

var
   rec_level       { recursion level }
      : integer;
   pvs_bug1, pvs_bug2, pvs_bug3 { bug in Pascal/VS argument passing }
      : real;

function simpson
         (function f(x : real) : real;
          a, da, fa, fm, fb, absarea, est, eps : real
         ) : real;

   {  This is the recursive Simpson's rule integration itself }

var
   dx, x1, x2, est1, est2, est3, f1, f2, f3, f4, sum,
   diff, tol, eps17, int1, int2, int3
      : real;

begin {------------------------    simpson    -----------------------}
   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);
   est1 := (fa + f1 + f2) * dx / 6.0;
   est2 := (f2 + fm + f3) * dx / 6.0;
   est3 := (f3 + f4 + fb) * dx / 6.0;
   absarea := absarea - abs(est) + abs(est1) + abs(est2) + abs(est3);
   sum := est1 + est2 + est3;
   rec_level := rec_level + 1;
   diff := abs(est - sum);
   tol := eps * absarea;
   if ((diff <=  tol) and (est <>  1.0))
         or (rec_level >=  MAXLEVEL) then
   begin
      simpson := sum;
      if rec_level >= MAXLEVEL then
         int_fail := int_fail + 1;
   end
   else
   begin
      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);
      simpson := int1 + int2 + int3;
   end;
   rec_level := rec_level - 1;
end { end simpson } ;

begin {-------------------------     integral    --------------------}
   rec_level := 1;
   int_fail := 0;
   pvs_bug1 := f(a);
   pvs_bug2 := 4.0 * f((a + b) / 2.0);
   pvs_bug3 := f(b);
   integral := simpson(f, a, b - a, pvs_bug1, pvs_bug2, pvs_bug3,
                       1.0, 1.0, eps);
end { end integral } ;