program integral1
   implicit none
   double precision x1,x2,s1,s2
   integer n,j
   x1=0.0d0
   x2=4.0d0*datan(1.0d0)
   n=8
   do j=1,7
      n=n*2
      call trapezoid(x1,x2,n,s1)
      call simpson(x1,x2,n,s2)
      write(*,100) n,s1,s2
   end do
   100 format(1x,'n=',1x,i4,1x,'trapezoidal',1x,f16.13,1x,'simpson',1x,f16.13)
   stop
   end

subroutine trapezoid(x1,x2,n,s)
   implicit none
   double precision x1,x2,s
   integer n
   double precision f,x,fx
   integer i
   s=0.0d0
   s=s+f(x1)
   do i=1,n-1
      x=x1+(x2-x1)/dble(n)*dble(i)
      fx=f(x)
      s=s+2.0d0*fx
   end do
   s=s+f(x2)
   s=s/2.0d0*(x2-x1)/dble(n)
   return
   end

subroutine simpson(x1,x2,n,s)
   implicit none
   double precision x1,x2,s
   integer n
   double precision f,xx1,xx2,fx1,fx2
   integer i
   s=0.0d0
   s=s+f(x1)
   do i=1,n/2
      xx1=x1+(x2-x1)/dble(n)*dble(i*2-1)
      fx1=f(xx1)
      xx2=x1+(x2-x1)/dble(n)*dble(i*2)
      fx2=f(xx2)
      s=s+4.0d0*fx1+2.0d0*fx2
   end do
   s=s-f(x2)
   s=s/3.0d0*(x2-x1)/dble(n)
   return
   end

function f(x)
   implicit none
   double precision f,x
   f=dsin(x)
   return
   end