C 台形公式を Richardson 加速してみる
      PROGRAM KASOKU
      DOUBLE PRECISION RMBERG,F,EPS,H,S,A,B
      EXTERNAL F
      EPS=0.5D-15
      A=0.0D0
      B=1.0D0
      H=0.1D0
      N=10
      WRITE(*,*)
     +  'Romberg 積分法による [0,1] 上の 1/(1+x) の定積分の計算'
      S=RMBERG(F,A,B,N,H,EPS)
      WRITE(*,200)'Romberg 加速法による値 :   ',S
      WRITE(*,*)'反復回数 : ',N,'最後の H : ',H
      WRITE(*,200)'真の値     :                  ',DLOG(2.0D0)
  200 FORMAT(1H ,A33,F18.15)
      END
C Romberg 法のサブルーチン
      DOUBLE PRECISION FUNCTION RMBERG(F,A,B,N,H,EPS)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION ROMB(100)
      EXTERNAL F
      IA=4
      I=1
   50 ROMB(I)=DAIKEI(F,A,B,N,H)
      M=IA
      DO 100 J=I-1,1,-1
       ROMB(J)=(M*ROMB(J+1)-ROMB(J))/(M-1)
       M=M*IA
  100 CONTINUE
      WRITE(*,200)(ROMB(J),J=1,I)
      IF ((I>1).AND.(DABS(ROMB(1)-ROMB(2))).LT.EPS) GO TO 170
      H=H/2
      N=N*2
      I=I+1
      GO TO 50
  170 RMBERG=ROMB(1)
      N=I
  200 FORMAT(1H ,8F18.15)
      END
C 台形公式のサブルーチン
      FUNCTION DAIKEI(F,A,B,N,H)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL F
      X=A
      DAIKEI=F(X)/2
      DO 100 I=1,N-1
       X=X+H
       DAIKEI=DAIKEI+F(X)
  100 CONTINUE
      DAIKEI=(DAIKEI+F(B)/2)*H
      END
C
      DOUBLE PRECISION FUNCTION F(X)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      F=1.0D0/(X+1)
      END