program ex0213b
      implicit real*8 (a-h,o-z)
      external f3
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      write(*,*) 'm = ? (n_max = 2**m)'
c      read(*,*)   m
c      write(*,*) 'm = ',m
c      write(*,*) 'n_max = ',2**m
      m=10
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
         call DEint(f3,a,b,n,sum)
         err = abs(sum-exa)
         write( *,800) j, dx, sum, exa, err
         write(16,800) j, dx, sum, exa, err
      end do			! 外側の DO LOOP の終り
c
 800   format(i4,5(1x,1pe12.5))
c
      stop
      end
c **********************************************************************
      subroutine DEint(FUNC,a,b,n,sum)
c **********************************************************************
c
c     This subroutine calculate the integral
c       sum = integral_a^b FUNC(x) dx
c     by using the Double Exponential formula (Takahashi-Mori).
c
c       sum= integral_{-infinity}^infinity FUNC(phi(t)) d(phi)/dt dt
c       phi = (b-a)/2*tanh(pi/2*sinh(t))+(b+a)/2
c
c **********************************************************************
      implicit real*8 (a-h,o-z)
      external FUNC
c
      tmax=4.7d0
      dt=tmax/n
      pi=4.0d0*atan(1.0d0)
c
      sum = 0.0d0
      do i = -n,n
         t=i*dt
         s=sinh(t)
         phi=(b-a)/2*tanh(pi/2*s)+(b+a)/2
         dphidt=pi/2*cosh(t)/(cosh(pi/2*s))**2*(b-a)/2
         sum=sum+FUNC(phi)*dphidt
      end do
         sum = sum*dt
C
      return
      end
c **********************************************************************
      function f3(x)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      eps=1.0d-16
c
      if(x.gt.eps) then
         f3 = -log(x)
      else
         f3 = -log(eps)
      endif
c
      end