C e^{-x^2} の定積分の近似計算
C EXTERNAL 文は函数を実引数とすることを可能にするための宣言文
      PROGRAM NORMAL
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(A=0.0D0,B=6.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)'台形公式による値     : ',S
      S=DAIKEI(G,A,B,N,H,1)
      WRITE(*,100)'同・逆からの和       : ',S
      S=SMPSON(G,A,B,N,H,0)
      WRITE(*,100)'SIMPSON 公式による値 : ',S
      S=SMPSON(G,A,B,N,H,1)
      WRITE(*,100)'同・逆からの和       : ',S
      WRITE(*,100)'真の値 (sqrt(π)/2)  :',DSQRT(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=DEXP(-X*X)
      RETURN
      END