PROGRAM KINKAI
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION NIBUN,NEWTON,LINITR
      EXTERNAL F,G
      S=NIBUN(F,1.0D0,2.0D0,I)
      WRITE(*,300) 'Bisection Method : ',S,'at I = ',I      
      S=NEWTON(F,G,1.0D0,2.0D0,I)
      WRITE(*,300) 'Newton'//CHAR(39)//'s Method  : ',S,'at I = ',I      
      S=LINITR(F,1.0D0,2.0D0,I)
      WRITE(*,300) 'Linear Iteration : ',S,'at I = ',I      
  300 FORMAT(1H ,A20,F18.15,A10,I7)
      END
C 
      DOUBLE PRECISION FUNCTION F(X)
      DOUBLE PRECISION X
      F=DEXP(X)-2*X-1
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION G(X)
      DOUBLE PRECISION X
      G=DEXP(X)-2
      RETURN
      END
C
C 二分法のサブルーチン
      DOUBLE PRECISION FUNCTION NIBUN(F,A,B,I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL F
      EPS=1.0D-16
      FA=F(A)
      FB=F(B)
      IF (FA*FB.GT.0.0D0) THEN
          WRITE(*,*) 'WRONG INTERVAL'
          RETURN
      END IF      
      XA=A
      XB=B
      DO 200 I=1,100  ! この反復回数は2^(-54)≒=10^(-16)より54以下で十分
          X=(XA+XB)/2
          WRITE(*,*)XA,X,XB,XB-XA
          FX=F(X)
          IF (FX*FA.LE.0.0D0) THEN
              XB=X
              FB=FX
          ELSE 
              XA=X
              FA=FX
          END IF
          IF (XB-XA.LT.EPS) GOTO 300
  200 CONTINUE 
  300 NIBUN=X
      RETURN
      END
C
C ニュートン法のサブルーチン
      DOUBLE PRECISION FUNCTION NEWTON(F,G,A,B,I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL F,G       ! G(X) = F'(X)
      EPS=1.0D-16
      NEWTON=B           ! 近似解の初期値
C     NEWTON=(A+B)/2     
      DO 100 I=1,100
          Y=F(NEWTON)/G(NEWTON)
          NEWTON=NEWTON-Y
          IF (DABS(Y).LT.EPS) RETURN
  100 CONTINUE
      RETURN
      END
C
C 線型反復法のサブルーチン
      DOUBLE PRECISION FUNCTION LINITR(F,A,B,I)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL F
      EPS=1.0D-16
      LINITR=B               ! 近似解の初期値
C      T=DEXP(B)-2            ! 使用する傾き
      T=(F(B)-F(A))/(B-A)
      DO 100 I=1,500
          Y=F(LINITR)/T
          IF (DABS(Y).LT.EPS) RETURN
          LINITR=LINITR-Y
  100 CONTINUE
      RETURN
      END