begin;
    comment  HAVIE INTEGRATOR.
    ALGORITHM 257, Robert N. Kubik,
    CACM 8 (1965) 381
    ;
    real A, B, EPS, MASK, Y, ANSWER;
    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;
                HAVIEINTEGRATOR := MASK;
                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;
    MASK := 9.99;
    ANSWER := HAVIEINTEGRATOR(Y, A, B, EPS, COS(Y), 12);
    NLCR;
    PRINTTEXT("integral(0,pi/2,cos(x)) = ");
    FLOT(10, 2, ANSWER);
    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)) = ");
    FLOT(10, 2, ANSWER);
    A := 1.0;
    B := 10.0;
    ANSWER := HAVIEINTEGRATOR(Y, A, B, EPS, LN(Y), 12);
    NLCR;
    PRINTTEXT("integral(1,10,ln(x)) = ");
    FLOT(10, 2, ANSWER);
    A := 0.0;
    B := 20.0;
    ANSWER := HAVIEINTEGRATOR(Y, A, B, EPS, Y  (1 ÷ 2) ÷ (EXP(Y - 4) + 1), 20);
    NLCR;
    PRINTTEXT("integral(1,20,x^(1/2)/(exp(x-4)+1)) = ");
    FLOT(10, 2, ANSWER);
    ;
end;