!###########################################################
! page.71 表4.1(例題1)
! 2分法のプログラム
!
! f(x)=x**3+6*x**2+21*x+32
!###########################################################
PROGRAM prog13
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: KMAX=99
!-------------
a=-3.d0
b=0.d0
c=0.d0
delta=1.d-10
eps=0.5d-5
!-------------
WRITE(*,"('a(0)=',1PE15.8,', b(0)=',1PE15.8,', delta=',1PE15.8,', eps=',1PE15.8)") a,b,delta,eps

fa=func(a)
fb=func(b)
IF(fa*fb>=0.d0)THEN
WRITE(*,*) 'initial value is incorrect' ; STOP
END IF

DO k=0,KMAX
    diff=dabs(a-b)
    c=(a+b)/2.d0
    fc=func(c)
    WRITE(*,"('k=',i2,', |a-b|=',1PE15.8,', c=',1PE15.8,', f(c)=',1PE15.8)") k,diff,c,fc

    IF(dabs(fc)<delta.or.diff<eps)GOTO 100

    IF(fc*fa>0.d0)THEN
    a=c ; fa=fc
    ELSE
    b=c ; fb=fc
    END IF
    
END DO
WRITE(*,*) '解が求まりませんでした'
STOP
100 CONTINUE
WRITE(*,"('解は x=',1PE15.8)") c
STOP
!###########################################################
CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION func(xx)
   REAL(8),INTENT(IN) :: xx
   REAL(8) :: func
   func=xx**3+6.d0*xx**2+21.d0*xx+32.d0
RETURN
END FUNCTION func
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog13