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