procedure kuncir(a, b, f, I, i, eps, N);
value a, b, eps, N;
integer N;
real a, b, I, i, eps;
real procedure f;
comment This procedure integrates the function f(x) using a
modified Simpson's Rule quadrature formula. The quadrature is
performed over j subintervals of [a,b] forming the total area I.
Convergence in each subintervals of length (b-a)2^n is indicated
when the relative difference between successive three-point and
five-point area approximations
A[3, j] = (b - a)*(g0 + 4*g2 + g4)/(3 * 2^(n + 1))
A[5, j] = (b - a)*(g0 + 4*g1 + 2*g2 + 4*g3 + g4)/(3 * 2^(n + 2))
is less than or equal to an appropriate portion of the over-all
tolerance eps (i.e., |A[5, j] - A[3, j]/A[5, j]| <= eps/2^n with n <= N.
`kuncir' will reduce the size of each interval until this
condition is satisfied.
Complete integration over [a,b] is indicated by i = b. A value
a <= i < b is indicates that the integration was terminated, leaving
I the true area under f in [a,b]. Further integration over [i,b]
will necessitate either the assignment of a larger N, a larger eps,
or an integral substitution reducing the slope of the integrand in
that interval. It is recommended that this procedure be used
between known integrand maxima and minima.;
begin integer m, n;
real d, h;
array g[0:4], A[0:2], S[1:N, 1:3];
I := 0;
i := 0;
m := 0;
n := 0;
g[0] := f(a);
g[2] := f((a + b)/2);
g[4] := f(b);
A[0] := (b - a)*(g[0] + 4*g[2] + g[4])/2;
AA: d := 2^n;
h := (b - a)/4/d;
g[1] := f(a + h*(4*m + 1));
g[3] := f(a + h*(4*m + 3));
A[1] := h*(g[0] + 4*g[1] + g[2]);
A[2] := h*(g[2] + 4*g[3] + g[4]);
if abs(((A[1] + A[2]) - A[0])/(A[1] + A[2])) > eps/d then 
begin 
  m := 2*m;
  n := n + 1;
  if n > N then goto CC;
  A[0] := A[1];
  S[n,1] := A[2];
  S[n,2] := g[3];
  S[n,3] := g[4];
  g[4] := g[2];
  g[2] := g[1];
  goto AA
end else 
begin 
  I := I + (A[1] + A[2])/3;
  m := m + 1;
  i := a + m*(b - a)/d;
  BB: if m = 2*(m % 2) then 
  begin 
    m := m % 2;
    n := n - 1;
    goto BB
  end;
  if (m != 1) | (n != 0) then 
  begin 
    A[0] := S[n,1];
    g[0] := g[4];
    g[2] := S[n,2];
    g[4] := S[n,3];
    goto AA
  end 
end;
CC: end;