program ex0213
      implicit real*8 (a-h,o-z)
c
c     This program calculate the integral
c       sum = integral_0^1 (-log(x))
c     and print out the numerical error
c     as a function of dividing number
c
c     NOTE: Since we cannot calculate log(0),
c           f(x) in this program returns
c              f(x) = - log(x)    (x >  eps)
c                   = - log(eps)) (x <= eps)
c           with eps=1.0d-16 (= 1.0 x 10^{-16})
c
      write(*,*) 'm = ? (n_max = 2**m)'
      read(*,*)   m
      write(*,*) 'm = ',m
      write(*,*) 'n_max = ',2**m
c      write(*,*) 'a, b = ?'
c      read(*,*)   a, b
c      write(*,*) 'a = ',a,' b = ',b
      a=0.0d0
      b=1.0d0
c
c Exact Value
      exa=1.0d0
c
      do j=1, m				! 外側の DO LOOP
             n=2**j
             sum1 = 0.0d0
             sum2 = 0.0d0
             sum3 = 0.0d0
             dx = (b - a)/n
          do i = 1, n			! 内側の DO LOOP
             x1 = a + (i - 1.0d0)*dx
             x2 = a + (i - 0.5d0)*dx
             x3 = a +  i         *dx
             sum1 = sum1 + f(x2)			! 中点
             sum2 = sum2 + (f(x1) + f(x3))/2		! 台形
             sum3 = sum3 + (f(x1) + 4*f(x2) + f(x3))/6	! シンプソン
          end do			! 内側の DO LOOP の終り
             sum1 = sum1*dx
             sum2 = sum2*dx
             sum3 = sum3*dx
c
             err1 = abs(sum1-exa)
             err2 = abs(sum2-exa)
             err3 = abs(sum3-exa)
c
         write( *,800) j, dx, exa, err1, err2, err3
         write(16,800) j, dx, exa, err1, err2, err3
c
      end do			! 外側の DO LOOP の終り
c
 800   format(i4,5(1x,1pe12.5))
c
      stop
      end
c **********************************************************************
      function f(x)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      eps=1.0d-16
c
      if(x.gt.eps) then
         f = -log(x)
      else
         f = -log(eps)
      endif
c
      end