Fix for ITP method

Fixes #25
readme update

README.md

@@ -58,28 +58,28 @@ To generate the documentation using [ford](https://github.com/Fortran-FOSS-Progr
 
 ### Methods
 
-The module contains the following methods (in alphabetical order):
+The module contains the following methods:
 
 Procedure | Description | Reference
 --- | --- | ---
-[`anderson_bjorck`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck.html) | Anderson-Bjorck method | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
-[`anderson_bjorck_king`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck_king.html) | a [variant](https://link.springer.com/content/pdf/bbm%3A978-3-642-05175-3%2F1.pdf) of `anderson_bjorck` | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
-[`barycentric`](https://jacobwilliams.github.io/roots-fortran/proc/barycentric.html) | Barycentric interpolation method | [Mendez & Castillo (2021)](https://www.researchgate.net/publication/352162661_A_highly_efficient_numerical_method_to_solve_non-linear_functions_using_barycentric_interpolation)
-[`bdqrf`](https://jacobwilliams.github.io/roots-fortran/proc/bdqrf.html) | Bisected Direct Quadratic Regula Falsi | [Gottlieb & Thompson (2010)](http://www.m-hikari.com/ams/ams-2010/ams-13-16-2010/gottliebAMS13-16-2010.pdf)
 [`bisection`](https://jacobwilliams.github.io/roots-fortran/proc/bisection.html) | Classic bisection method | Bolzano (1817)
-[`blendtf`](https://jacobwilliams.github.io/roots-fortran/proc/blendtf.html) | Blended method of trisection and false position | [Badr, Almotairi, Ghamry (2021)](https://www.mdpi.com/2227-7390/9/11/1306/htm)
+[`regula_falsi`](https://jacobwilliams.github.io/roots-fortran/proc/regula_falsi.html) | Classic regula falsi method | ?
+[`muller`](https://jacobwilliams.github.io/roots-fortran/proc/muller.html) | Improved Muller method (for real roots only) | [Muller (1956)](https://www.ams.org/journals/mcom/1956-10-056/S0025-5718-1956-0083822-0/S0025-5718-1956-0083822-0.pdf)
 [`brent`](https://jacobwilliams.github.io/roots-fortran/proc/brent.html) | Classic Brent's method (a.k.a. Zeroin) | [Brent (1971)](https://maths-people.anu.edu.au/~brent/pd/rpb005.pdf)
 [`brenth`](https://jacobwilliams.github.io/roots-fortran/proc/brenth.html) | SciPy variant of `brent` | [SciPy](https://github.com/scipy/scipy/)
 [`brentq`](https://jacobwilliams.github.io/roots-fortran/proc/brentq.html) | SciPy variant of `brent` | [SciPy](https://github.com/scipy/scipy/)
-[`chandrupatla`](https://jacobwilliams.github.io/roots-fortran/proc/chandrupatla.html) | Hybrid quadratic/bisection algorithm | [Chandrupatla (1997)](https://dl.acm.org/doi/10.1016/S0965-9978%2896%2900051-8)
 [`illinois`](https://jacobwilliams.github.io/roots-fortran/proc/illinois.html) | Illinois method | [Dowell & Jarratt (1971)](https://personal.math.ubc.ca/~loew/mech2/Dowell+Jarratt.pdf)
-[`itp`](https://jacobwilliams.github.io/roots-fortran/proc/itp.html) | Interpolate Truncate and Project method | [Oliveira & Takahashi (2020)](https://dl.acm.org/doi/abs/10.1145/3423597)
-[`muller`](https://jacobwilliams.github.io/roots-fortran/proc/muller.html) | Improved Muller method (for real roots only) | [Muller (1956)](https://www.ams.org/journals/mcom/1956-10-056/S0025-5718-1956-0083822-0/S0025-5718-1956-0083822-0.pdf)
 [`pegasus`](https://jacobwilliams.github.io/roots-fortran/proc/pegasus.html) | Pegasus method | [Dowell & Jarratt (1972)](https://link.springer.com/article/10.1007/BF01932959)
-[`regula_falsi`](https://jacobwilliams.github.io/roots-fortran/proc/regula_falsi.html) | Classic regula falsi method | ?
+[`anderson_bjorck`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck.html) | Anderson-Bjorck method | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
+[`anderson_bjorck_king`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck_king.html) | a [variant](https://link.springer.com/content/pdf/bbm%3A978-3-642-05175-3%2F1.pdf) of `anderson_bjorck` | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
 [`ridders`](https://jacobwilliams.github.io/roots-fortran/proc/ridders.html) | Classic Ridders method | [Ridders (1979)](https://cs.fit.edu/~dmitra/SciComp/Resources/RidderMethod.pdf)
 [`toms748`](https://jacobwilliams.github.io/roots-fortran/proc/toms748.html) | Algorithm 748 | [Alefeld, Potra, Shi (1995)](https://dl.acm.org/doi/abs/10.1145/210089.210111)
+[`chandrupatla`](https://jacobwilliams.github.io/roots-fortran/proc/chandrupatla.html) | Hybrid quadratic/bisection algorithm | [Chandrupatla (1997)](https://dl.acm.org/doi/10.1016/S0965-9978%2896%2900051-8)
+[`bdqrf`](https://jacobwilliams.github.io/roots-fortran/proc/bdqrf.html) | Bisected Direct Quadratic Regula Falsi | [Gottlieb & Thompson (2010)](http://www.m-hikari.com/ams/ams-2010/ams-13-16-2010/gottliebAMS13-16-2010.pdf)
 [`zhang`](https://jacobwilliams.github.io/roots-fortran/proc/zhang.html) | Zhang's method (with corrections from Stage) | [Zhang (2011)](https://www.cscjournals.org/download/issuearchive/IJEA/Volume2/IJEA_V2_I1.pdf)
+[`itp`](https://jacobwilliams.github.io/roots-fortran/proc/itp.html) | Interpolate Truncate and Project method | [Oliveira & Takahashi (2020)](https://dl.acm.org/doi/abs/10.1145/3423597)
+[`barycentric`](https://jacobwilliams.github.io/roots-fortran/proc/barycentric.html) | Barycentric interpolation method | [Mendez & Castillo (2021)](https://www.researchgate.net/publication/352162661_A_highly_efficient_numerical_method_to_solve_non-linear_functions_using_barycentric_interpolation)
+[`blendtf`](https://jacobwilliams.github.io/roots-fortran/proc/blendtf.html) | Blended method of trisection and false position | [Badr, Almotairi, Ghamry (2021)](https://www.mdpi.com/2227-7390/9/11/1306/htm)
 
 In general, all the methods are guaranteed to converge. Some will be more efficient (in terms of number of function evaluations) than others for various problems. The methods can be broadly classified into three groups:
 

src/root_module.F90

@@ -1237,7 +1237,7 @@ end subroutine bdqrf
 ! * D. E. Muller, "A Method for Solving Algebraic Equations Using an Automatic Computer",
 ! Mathematical Tables and Other Aids to Computation, 10 (1956), 208-215.
 ! [link](https://www.ams.org/journals/mcom/1956-10-056/S0025-5718-1956-0083822-0/S0025-5718-1956-0083822-0.pdf)
-! * [Roots.jl](https://github.com/JuliaMath/Roots.jl/blob/master/src/simple.jl), 
+! * [Roots.jl](https://github.com/JuliaMath/Roots.jl/blob/master/src/simple.jl),
 ! Julia version of standard Muller
 
  subroutine muller (me,ax,bx,fax,fbx,xzero,fzero,iflag)
@@ -2517,24 +2517,25 @@ subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)
  integer :: j !! iteration counter
  logical :: root_found !! convergence in x
  logical :: fail !! if we can't interpolate
+ real(wp) :: factor !! factor to ensure f(ax)<f(bx)
 
  real(wp),parameter :: log2 = log(2.0_wp)
 
  ! initialize:
- if (fax < fbx) then
- a = ax
- b = bx
- ya = fax
- yb = fbx
+ if (fax < fbx) then ! a requirement for the algorithm ?
+ factor = 1.0_wp
  else
- a = bx
- b = ax
- ya = fbx
- yb = fax
+ ! flip the sign of the function to ensure fax < fbx
+ ! note that a < b is already enforced by the calling routine.
+ factor = -1.0_wp
  end if
+ a = ax
+ b = bx
+ ya = factor*fax
+ yb = factor*fbx
  iflag = 0
  term = (b-a)/(2.0_wp*me%rtol)
- n12 = ceiling ( log(term) / log2 ) ! ceiling(log2(term))
+ n12 = ceiling ( log(abs(term)) / log2 ) ! ceiling(log2(term))
  nmax = n12 + me%n0
  aprev = huge(1.0_wp) ! initialize to unusual values
  bprev = huge(1.0_wp) !
@@ -2571,7 +2572,7 @@ subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)
  end if
 
  ! updating interval:
- yitp = me%f(xitp)
+ yitp = factor*me%f(xitp)
  if (me%solution(xitp,yitp,xzero,fzero)) return
  if (yitp > 0.0_wp) then
  b = xitp
@@ -2591,7 +2592,7 @@ subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)
  ! [this can happen in the test cases].
  ! So just do a bisection
  xitp = bisect(a,b)
- yitp = me%f(xitp)
+ yitp = factor*me%f(xitp)
  if (me%solution(xitp,yitp,xzero,fzero)) return
  if (ya*yitp<0.0_wp) then
  ! root lies between a and xitp

test/root_tests.f90

@@ -1300,12 +1300,27 @@ subroutine problems(x, ax, bx, fx, xroot, cases, num_of_problems, latex)
  if (present(x)) f = x
  if (present(latex)) latex = 'x'
 
+ ! another linear test case (see Issue #25)
+ case(133)
+
+ tmp : block
+ real(wp),parameter :: x1 = 0.0_wp
+ real(wp),parameter :: x2 = 157.08_wp
+ real(wp),parameter :: y1 = 1012.0_wp
+ real(wp),parameter :: y2 = -1114.0_wp
+ a = x1
+ b = x2
+ root = 7.4771853245531515E+01_wp
+ if (present(x)) f = y1 + ((x - x1) / (x2 - x1)) * (y2 - y1)
+ if (present(latex)) latex = '1012 + (x/x2)(-2126)'
+ end block tmp
+
  case default
  write(*,*) 'invalid case: ', nprob
  error stop 'invalid case'
  end select
 
- if (present(num_of_problems)) num_of_problems = 132
+ if (present(num_of_problems)) num_of_problems = 133
 
  ! outputs:
  if (present(ax)) ax = a