program trap3
c
      integer out
      real sum, upper, lower, tol
      external f, df
      common /inout/ out
      data lower/1.0/, upper/9.0/, tol/1.0e-5/
c
      out = 6
      write(out, 101)
      call trapez(lower, upper, tol, sum, f, df)
      write(out, 104) sum
      stop
101   format(/' trapezoidal integration with end correction')
104   format(/' area =', f10.5/)
      end
      function f(x)
c
      f = 1.0 / x
      return
      end
      function df(x)
c
      df = -1.0/(x * x)
      return
      end
      subroutine trapez(lower, upper, tol, sum, f, df)
c
c -- numerical integration by the trapezoidal method
c
      integer out, pieces, i, p2
      real x, delta, lower, upper, sum, tol
      real endsum, midsum, sum1, endcor
      common /inout/ out
c
      pieces  = 1
      delta = (upper - lower) / pieces
      endsum = f(lower) + f(upper)
      endcor = (df(upper) - df(lower)) / 12.0
      sum = endsum * delta / 2.0
      write(out, 101) sum
      midsum = 0.0
5     pieces = pieces * 2
       p2 = pieces / 2
       sum1 = sum
       delta = (upper - lower) / pieces
       do 10 i = 1, p2
        x = lower + delta *(2 * i - 1)
        midsum = midsum + f(x)
10     continue
       sum = (endsum + 2.0*midsum) * delta * 0.5
     *   - delta * delta * endcor
       write(out, 102) pieces, sum
      if (abs(sum - sum1) .gt. abs(tol * sum)) goto 5
      return
101   format(/'       1', f9.5)
102   format(1x, i7, f9.5)
      end