C 前進差分公式を Richardson 加速してみる
      PROGRAM KASOKU
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION B(20)
      F(X)=DSIN(X)
      DF(X,H)=(F(X+H)-F(X))/H
C
      WRITE(*,*)'sin(x) の x=1 における微分を'
      WRITE(*,*)'前進差分と Richardson 加速で計算します．'
      X=1.0D0
C
      EPS=0.00000000000001
      IA=2
      H=0.1D0
      N=10
      DO I=1,N
        B(I)=DF(X,H)
        M=IA
        DO J=I-1,1,-1
          B(J)=(M*B(J+1)-B(J))/(M-1)
          M=M*IA
        END DO
        WRITE(*,200)(B(J),J=1,I)
        IF ((I>1).AND.(DABS(B(1)-B(2))).LT.EPS) GO TO 170
        H=H/2
      END DO
  170 WRITE(*,300)'Richardson 加速法による値 :   ',B(1)
C     WRITE(*,*)'停止ε : ',EPS,'反復回数 : ',I,'最後の H : ',H
      WRITE(*,*)'Stopping Eps : ',EPS,', Rep. No : ',I,', Last H : ',H
      WRITE(*,300)'真の値     :                  ',DCOS(1.0D0)
  200 FORMAT(1H ,8F18.15)
  300 FORMAT(1H ,A33,F18.15)
      END