program prog023
      implicit real*8 (a-h,o-z)
c
c     This program calculate the integral
c       sum = integral_a^b 1/sqrt(x) dx
c     and print out the numerical error
c     as a function of dividing number
c
      write(*,*) 'm = ? (n_max = 2**m)'
      read(*,*)   m
      write(*,*) 'm = ',m
      write(*,*) 'n_max = ',2**m
      write(*,*) 'a, b = ?'
      read(*,*)   a, b
      write(*,*) 'a = ',a,' b = ',b
c
c Exact Value
      exa=2*(sqrt(b)-sqrt(a))
c
      write(*,*) 'j, dx, Exact, Err1, Err2, Err3'	! 何を書き出すのか、
      write(16,*) 'j, dx, Exact, Err1, Err2, Err3'	! メモしておく。
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)
      f = 1.0d0/sqrt(x)
      end