C int 1/(x^2+1)dx calculated by exponential transform.
      PROGRAM INFINT
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(A=0.0D0,B=20.0D0,PI=3.14159265358979323846D0)
      EXTERNAL G
   10 WRITE(*,*)'分割数は？（偶数で） : '
      READ(*,*)N
      H=(B-A)/N
      WRITE(*,100)'刻み幅 =               ',H
      S=DAIKEI(G,A,B,N,H,0)
      WRITE(*,100)'台形公式による値     : ',2*S
      S=DAIKEI(G,A,B,N,H,1)
      WRITE(*,100)'同・逆からの和       : ',2*S
      S=SMPSON(G,A,B,N,H,0)
      WRITE(*,100)'SIMPSON 公式による値 : ',2*S
      S=SMPSON(G,A,B,N,H,1)
      WRITE(*,100)'同・逆からの和       : ',2*S
      WRITE(*,100)'真の値 (π/2)  :',PI/2
      GO TO 10
  100 FORMAT(1H ,A23,D22.15)
      END
C 台形公式のサブルーチン
      FUNCTION DAIKEI(F,A,B,N,H,IND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL F
      IF (IND.EQ.0) THEN
          X=A
          HH=H
      ELSE
          X=B
          HH=-H
      END IF
      DAIKEI=F(X)/2
      DO 100 I=1,N-1
          X=X+HH
          DAIKEI=DAIKEI+F(X)
  100 CONTINUE
          X=X+HH
      DAIKEI=(DAIKEI+F(X)/2)*H
      RETURN
      END
C SIMPSON 公式のサブルーチン
      FUNCTION SMPSON(F,A,B,N,H,IND)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL F
      IF (IND.EQ.0) THEN
          X=A
          HH=H
      ELSE
          X=B
          HH=-H
      END IF
      SMPSON=F(X)
      DO 100 I=1,N-2,2
       X=X+HH
       SMPSON=SMPSON+F(X)*4
       X=X+HH
       SMPSON=SMPSON+F(X)*2
  100 CONTINUE
      X=X+HH
      SMPSON=SMPSON+F(X)*4
      X=X+HH
      SMPSON=(SMPSON+F(X))*H/3
      RETURN
      END
C
C 被積分函数のサブルーチン
      FUNCTION G(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      G=1.0D0/(DEXP(-X)+DEXP(X))
      RETURN
      END