program simpson
   implicit none
   integer,parameter::n=20
   double precision,parameter::x1=-1.0d0
   double precision,parameter::x2=1.0d0
   double precision f,x,dx,s1
   double precision exact,epsilon
   integer i

   f(x)=sqrt(1.d0-x**2)

   dx=(x2-x1)/dble(n)
   s1=f(x1)+f(x2)
   do i=1,n-1,2
      x=dble(i)*dx-1.d0
      s1=s1+4.d0*f(x)
   end do
   do i=2,n-2,2
      x=dble(i)*dx-1.d0
      s1=s1+2.d0*f(x)
   end do
   s1=s1*dx/3.d0
   exact=atan(1.d0)*2.d0
   epsilon=abs(s1-exact)

   write(*,100) s1,epsilon,epsilon/exact
   100 format(1x,'s1= ',1pe14.7,' err= ',1pe14.7,' relerr= ', 1pe14.7)
   stop
   end