PROGRAM teste
c     floating point, double precision arithmetic test
C     MAIN PROGRAM - EVALUATES INTEGRAL TO AN ACCURACY OF 8D
C                    INCREASING THE NUMBER OF STRIPS UNTIL REQUIRED
C                    ACCURACY IS REACHED
C               (A)  USING TRAPEZOIDAL RULE
C               (B)  USING SIMPSON'S RULE
C
      implicit double precision (a-h,o-z)
      character*18 type
C     (A)
C
      type='Trapezoidal Rule  '
      A=1.0
      B=2.2
C
      N=2
      ICOUNT=0
      WRITE(unit=6,fmt=103)type
   10 W=(B-A)/N
      ICOUNT=ICOUNT+1
      AREA=TRAP(A,B,N,W)
      IF(ICOUNT.EQ.200) THEN
         ICOUNT=0
         WRITE(unit=6,fmt=102)N,AREA
      END IF
      IF(N.NE.2) THEN
         IDIFF=ABS(AREA-ALAST)*100000000.0
         IF(IDIFF.eq.0) THEN
            N2=N*2
            W2=(B-A)/N2
            AREA2=TRAP(A,B,N2,W2)
            ICHECK=ABS(AREA2-AREA)*100000000.0
            IF(ICHECK.eq.0) THEN
                 write(unit=6,fmt=101)a,b,area,n
                 go to 49
            end if
         end if
      end if
      N=N+1
      ALAST=AREA
      GO TO 10
C
  101 FORMAT(1H0,/' Integral of function between the limits ',
     1 F5.2,'  and  ',F5.2/10X,'=',F14.8/
     2 ' Required accuracy reached by using minimum of ',
     3 I5,' strips' )
  102 FORMAT(1H ,'No. of strips =  ',I5,10X,
     1 'Area under curve =  ',F14.9)
  103 FORMAT(1H1,'Numerical Integration Using -  ',
     1 a18,//)
C
   49 CONTINUE
C     (B)
      type='Simpson''s Rule  '
C
      N=2
      icount=0
      WRITE(unit=6,fmt=103)type
  210 W=(B-A)/N
      icount=icount+1
      AREA=SIMP(A,B,N,W)
      if(icount.eq.5) then
         icount=0
         WRITE(unit=6,fmt=102)N,AREA
      end if
      IF(N.ne.2) THEN
         IDIFF=ABS(AREA-ALAST)*100000000.0
         IF(IDIFF.eq.0) then
            N2=N*2
            W2=(B-A)/N2
            AREA2=SIMP(A,B,N2,W2)
            ICHECK=ABS(AREA2-AREA)*100000000.0
            IF(ICHECK.eq.0) then
               write(unit=6,fmt=101)a,b,area,n
               stop
            end if
         end if
      end if
      N=N+2
      ALAST=AREA
      GO TO 210
C
      END
C     INTEGRATION OF FUNCTION USING TRAPEZOIDAL RULE
C
      DOUBLE PRECISION FUNCTION TRAP(A,B,N,W)
      implicit double precision (a-h,o-z)
      F(X)=1+LOG(X)
      SUM=(F(A)+F(B))/2
      N1=N-1
      DO 2 I=1,N1
    2 SUM=SUM+F(A+W*I)
      TRAP=SUM*W
      RETURN
      END
C
C     INTEGRATION OF FUNCTION USING SIMPSON'S RULE
C
      DOUBLE PRECISION FUNCTION SIMP(A,B,N,W)
      implicit double precision (a-h,o-z)
      F(X)=1+LOG(X)
      SUM=F(A)+F(B)+4*F(B-W)
      IF(N.EQ.2) GO TO 3
      N3=N-3
      DO 2 I=1,N3,2
    2 SUM=SUM+4*F(A+W*I)+2*F(A+W*(I+1))
    3 SIMP=SUM*W/3.0
      RETURN
      END