!###########################################################
! page.71 表4.1(例題1)
! 2分法のプログラム(Fortran77)
!
! f(x)=x**3+6*x**2+21*x+32
!###########################################################
      program prog13
      parameter (kmax=99)
      implicit double precision(a-h,o-z)

      a=-3.d0
      b=0.d0
      c=0.d0
      delta=1.d-10
      eps=0.5d-5
      
      write(*,*) 'a(0)=',a,', b(0)=',b,', delta=',delta,'eps=',eps
      fa=func(a)
      fb=func(b)
       if(fa*fb>=0.d0) then
        write(*,*) 'initial value is incorrect' 
        stop
       end if

      do 100 k=0,kmax
       diff=dabs(a-b)
       c=(a+b)/2.d0
       fc=func(c)
       write(*,900) k,diff,c,fc
  900  format('k=',i4,' |a-b|=',1pe13.6,' c=',1pe13.6,' f(c)=',1pe13.6)
       if(dabs(fc)<delta.or.diff<eps) then 
            write(*,*) '解は x=',c,'です。'
            stop
       end if
       if(fc*fa>0.d0) then
           a=c ; fa=fc
         else
           b=c ; fb=fc
       end if
  100 continue
      write(*,*) '解が求まりませんでした。'
      stop
      end
      function func(xx)
       real*8 xx,func
       func=xx**3+6.d0*xx**2+21.d0*xx+32.d0
      return
      end