      program prog022
      implicit real*8 (a-h,o-z)
c
c     This program calculate the integral
c       sum = integral_a^b 1/sqrt(x) dx
c     by using the Simpson formula.
c
      write(*,*) 'n, a, b = ?'
      read(*,*)   n, a, b
      write(*,*) 'n = ', n,' a = ',a,' b = ',b
c
         sum = 0.0d0
         dx = (b - a)/n
c
      do i = 1, n
         x1 = a + (i - 1.0d0)*dx	! 短冊の左端
         x2 = a + (i - 0.5d0)*dx	! 短冊の中央
         x3 = a +  i         *dx	! 短冊の右端
         sum = sum + (f(x1) + 4*f(x2) + f(x3))/6	! 1:4:1 の比で足す。
      end do
         sum = sum*dx
c
      exa=2*(sqrt(b)-sqrt(a))
      write(*,*) 'Sum    = ',sum
      write(*,*) 'Exact  = ',exa
      stop
      end
c **********************************************************************
      function f(x)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      f = 1.0d0/sqrt(x)
      end