!=======================================================================
! ch11_false_position.f90
!
! Find one root of f(x) = 0 by the false position (regula falsi) method.
!
! Origin:
!   Based on the false-position routines by Alexander Godunov in
!   Root5.m and Roots5.cpp (February 2022). This Fortran program uses
!   the same ordinary false-position interpolation, with the layout
!   and comments of the Chapter 11 bisection program.
!
! Example:
!   f(x) = sqrt(x) - cos(x) on [0,1].
!   Reference root = 0.64171437087288... .
!
! The program reports the iteration count, root estimate, and f(root).
! The requested tolerance applies to the residual |f(root)|. Unlike
! bisection, ordinary false position may leave one endpoint fixed for
! many iterations, so the bracket width is not the stopping test.
!=======================================================================
program main
implicit none
double precision f,x1,x2,root,eps
integer flag
external f

! Initial bracket and requested residual tolerance
x1  = 0.0d0
x2  = 1.0d0
eps = 1.0d-6

write(*,100)
call false_position(f,x1,x2,eps,root,flag)

! Report a missing bracket or invalid input before printing a root
if(flag == 0) then
  write(*,104)
  stop
end if

write(*,101) abs(flag),eps
if(flag < 0) write(*,103)
write(*,102) root

! Check the residual of the reported estimate
write(*,105) f(root)

100 format(' A single root of f(x) ',/, &
           ' Method - False position')
101 format(' iterations        = ',i3,/,&
           ' residual tolerance = ',1pe12.5)
102 format(' root              = ',1pe12.5)
103 format(' iteration limit reached before the requested residual')
104 format(' invalid interval or no initial sign change')
105 format(' f(root)           = ',1pe12.5)

end program main


subroutine false_position(f,x1,x2,eps,Root,flag)
!=======================================================================
! Locate one root of f(x) on [x1,x2] by ordinary false position.
!
! Input:
!   f     - function to evaluate, continuous on the initial bracket
!   x1    - left endpoint of the initial bracket
!   x2    - right endpoint of the initial bracket
!   eps   - requested maximum residual |f(Root)| (positive)
!
! Output:
!   Root  - most recent false-position estimate when flag is nonzero
!   flag  - convergence indicator and iteration count
!           flag > 0: |f(Root)| <= eps; flag is the iteration count
!           flag = 0: invalid interval/tolerance or no sign change
!           flag < 0: iteration limit reached; abs(flag) is the count
!
! Notes:
!   A sign change guarantees a root only when f is continuous.
!   The maximum number of iterations is 200. One bracket endpoint
!   can remain nearly fixed, so a small interval width is not assumed.
!   An endpoint root is returned immediately with a positive flag.
!=======================================================================
implicit none
double precision f,x1,x2,eps,Root
double precision a,b,c,fa,fb,fc
integer i,flag
integer, parameter :: iter=200
external f

! Check that the inputs describe a usable initial bracket
flag = 0
if(x1 >= x2 .or. eps <= 0.0d0) return

a = x1
b = x2
fa = f(a)
fb = f(b)
if(fa*fb > 0.0d0) return

! Handle an exact zero at an endpoint before forming a secant line.
! Flag 1 follows the existing convention that positive means success.
if(fa == 0.0d0) then
  Root = a
  flag = 1
  return
end if
if(fb == 0.0d0) then
  Root = b
  flag = 1
  return
end if

! Interpolate between the current endpoints and retain the sign change.
do i=1,iter
  c = b - fb*(b-a)/(fb-fa)
  fc = f(c)

  ! False position uses the residual, not the bracket width, to stop.
  if(abs(fc) <= eps) then
    Root = c
    flag = i
    return
  end if

  if(fa*fc <= 0.0d0) then
    b = c
    fb = fc
  else
    a = c
    fa = fc
  end if
end do

Root = c
flag = -iter
end subroutine false_position


!-----------------------------------------------------------------------
function f(x)
implicit none
double precision f,x
! Example function. The chosen bracket must stay in the domain x >= 0.
! To try the alternative from the original programs, use [-1,0].
! f = x + cos(x)
f = sqrt(x) - cos(x)
end function f