!=======================================================================
! ch11_multiple_roots.f90
!
! Find several distinct roots by scanning an interval and using
! bisection on each subinterval with a sign change.
!
! Origin:
!   Based on BForce.f90 by Alexander Godunov (January 2010).
!   The original default was a brute-force scan followed by bisection.
!   That method is retained here without the unused method switch.
!
! Example:
!   f(x) = exp(x) - x*sin(pi*x/2) on [-10.5,2].
!   The chapter example has six distinct roots in this interval.
!
! The tolerance applies to the final width of each root bracket.
! A sign-change scan can miss even-multiplicity roots, pairs of roots
! within one subinterval, and roots near a discontinuity. Inspect a
! plot and repeat with more subintervals when counting roots matters.
!=======================================================================
program main
implicit none
integer, parameter :: n=1000
double precision f,x1,x2,eps,roots(n+1)
integer nroots,nfail,i
external f

! Search interval, number of scan subintervals, and bracket-width tolerance
x1  = -10.5d0
x2  =   2.0d0
eps =   1.0d-7

write(*,100) x1,x2,n,eps
call MultipleRoots(f,x1,x2,n,eps,roots,nroots,nfail)

if(nroots == 0) then
  write(*,101)
else
  write(*,102)
  do i=1,nroots
    write(*,103) i,roots(i),f(roots(i))
  end do
end if
if(nfail > 0) write(*,104) nfail

100 format(' Roots of f(x) by interval scanning and bisection',/,&
           ' interval             = [',f7.3,',',f7.3,']',/,&
           ' scan subintervals    = ',i5,/,&
           ' bracket-width tol.   = ',1pe12.5)
101 format(' no roots detected by the scan')
102 format(' number       root                  f(root)')
103 format(i5,3x,f18.10,3x,1pe16.8)
104 format(' brackets not refined to the requested width: ',i5)
end program main


subroutine MultipleRoots(f,x1,x2,n,eps,roots,nroots,nfail)
!=======================================================================
! Scan [x1,x2] and bisect subintervals with opposite endpoint signs.
!
! Input:
!   f       - function continuous on the search interval
!   x1,x2   - endpoints of the search interval, with x1 < x2
!   n       - number of equal scan subintervals (positive)
!   eps     - maximum final bracket width (positive)
!
! Output:
!   roots   - detected roots, in increasing order; capacity n+1
!   nroots  - number of roots stored
!   nfail   - number of sign-change brackets that reached the limit
!
! Notes:
!   A zero at a grid point is recorded once, including at x2.
!   A sign change only certifies a root when f is continuous there.
!   The scan does not certify that all roots have been found.
!=======================================================================
implicit none
integer n,nroots,nfail,j,flag
double precision f,x1,x2,eps,roots(n+1)
double precision a,b,fa,fb,dx,root
external f

nroots = 0
nfail = 0
if(n <= 0 .or. x1 >= x2 .or. eps <= 0.0d0) return

dx = (x2-x1)/dble(n)
a = x1
fa = f(a)

! Each grid point is checked once; a sign change brackets another root.
do j=1,n
  b = x1 + dble(j)*dx
  fb = f(b)

  if(fa == 0.0d0) then
    nroots = nroots+1
    roots(nroots) = a
  else if(fa*fb < 0.0d0) then
    call BisectBracket(f,a,b,eps,root,flag)
    if(flag > 0) then
      nroots = nroots+1
      roots(nroots) = root
    else
      nfail = nfail+1
    end if
  end if

  a = b
  fa = fb
end do

! The rightmost grid point is not a left endpoint in the loop.
if(fa == 0.0d0) then
  nroots = nroots+1
  roots(nroots) = a
end if
end subroutine MultipleRoots


subroutine BisectBracket(f,x1,x2,eps,root,flag)
!=======================================================================
! Refine one sign-change bracket by bisection.
!
! Input:
!   f       - function continuous on [x1,x2]
!   x1,x2   - endpoints with opposite signs of f
!   eps     - requested maximum final bracket width
!
! Output:
!   root    - midpoint of the final bracket
!   flag    - positive iteration count on convergence; -200 on limit
!
! Notes:
!   The stopping rule is bracket width, as in ch11_bisection.f90.
!=======================================================================
implicit none
integer flag,i
integer, parameter :: iter=200
double precision f,x1,x2,eps,root,a,b,c,fa,fc
external f

a = x1
b = x2
fa = f(a)

do i=1,iter
  c = a + (b-a)/2.0d0
  fc = f(c)
  if(fa*fc <= 0.0d0) then
    b = c
  else
    a = c
    fa = fc
  end if
  if(abs(b-a) <= eps) exit
end do

root = a + (b-a)/2.0d0
if(i > iter) then
  flag = -iter
else
  flag = i
end if
end subroutine BisectBracket


function f(x)
!-----------------------------------------------------------------------
! Example from Chapter 11. Replace this expression for another problem.
! Choose a scan interval on which the new function is continuous.
!-----------------------------------------------------------------------
implicit none
double precision f,x
double precision, parameter :: pi=3.141592653589793d0
f = exp(x) - x*sin(pi*x/2.0d0)
end function f