```begin;
comment  HAVIE INTEGRATOR.
ALGORITHM 257, Robert N. Kubik,
CACM 8 (1965) 381
;
integer COUNT;
real procedure HAVIEINTEGRATOR(X, A, B, EPS, INTEGRAND, M);
value A, B, EPS, M;
integer M;
real INTEGRAND, X, A, B, EPS;
comment  This algorithm performs numerical integration of
definite integrals using an equidistant sampling of the
function and repeated halving of the sampling interval.
Each halving allows the calculation of a trapezium and
a tangent formula on a finer grid, but also the calcul-
ation of several higher order formulas which are defined
implicitly. The two families of approximate solutions
will normally bracket the value of the integral and from
these convergence is tested on each of the several orders
of approximation. The algorithm is based on a private
communication from F. Haavie of the Institutt for Atom-
energi Kjeller Research Establishment, Norway. A Fortran
version is in use on the Philco-2000. ...;
begin;
real H, ENDPTS, SUMT, SUMU, D;
integer I, J, K, N;
real array T, U, TPREV, UPREV[1 : M];
X := A;
ENDPTS := INTEGRAND;
COUNT := 1;
X := B;
ENDPTS := 0.5  (INTEGRAND + ENDPTS);
COUNT := COUNT + 1;
SUMT := 0.0;
I := N := 1;
H := B - A;
ESTIMATE: T[1] := H  (ENDPTS + SUMT);
SUMU := 0.0;
comment  t[1] = h*(0.5*f[0]+f[1]+f[2]+...+0.5*f[2^(i-1)])
;
X := A - H ÷ 2.0;
for J := 1 step 1 until N do
begin;
X := X + H;
SUMU := SUMU + INTEGRAND;
COUNT := COUNT + 1;
end ;
U[1] := H  SUMU;
K := 1;
comment  u[1] = h*(f[1/2]+f[3/2]+...f[(2^i-1)/2],
k corresponds to approximate solution with truncation
error term of order 2k
;
TEST: if ABS(T[K] - U[K]) < EPS then begin;
HAVIEINTEGRATOR := 0.5  (T[K] + U[K]);
goto EXIT;
end ;
if K  I then begin;
D := 2  (2  K);
T[K + 1] := (D  T[K] - TPREV[K]) ÷ (D - 1.0);
TPREV[K] := T[K];
U[K + 1] := (D  U[K] - UPREV[K]) ÷ (D - 1.0);
UPREV[K] := U[K];
comment  This implicit formulation of the higher
order integration formulas is given in
[ROMBERG, W. ...
;
K := K + 1;
if K = M then begin;
goto EXIT;
end ;
goto TEST;
end ;
H := H ÷ 2.0;
SUMT := SUMT + SUMU;
TPREV[K] := T[K];
UPREV[K] := U[K];
I := I + 1;
N := 2  N;
goto ESTIMATE;
EXIT: NLCR;
NLCR;
PRINTTEXT("i: ");
ABSFIXT(3, 0, I);
PRINTTEXT("k: ");
ABSFIXT(3, 0, K);
PRINTTEXT(" count:");
ABSFIXT(7, 0, COUNT);
end HAVIEINTEGRATOR;
comment  Following is a driver program to test havieintegrator
;
A := 0.0;
B := 1.5707963;
EPS := 0.00001;
ANSWER := HAVIEINTEGRATOR(Y, A, B, EPS, COS(Y), 12);
NLCR;
PRINTTEXT("integral(0,pi/2,cos(x)) = ");
A := 0.0;
B := 4.3;
ANSWER := HAVIEINTEGRATOR(Y, A, B, EPS, EXP(-Y  Y), 12);
NLCR;
PRINTTEXT("integral(0,4.3,exp(-x*x)) = ");
A := 1.0;
B := 10.0;
ANSWER := HAVIEINTEGRATOR(Y, A, B, EPS, LN(Y), 12);
NLCR;
PRINTTEXT("integral(1,10,ln(x)) = ");