Handle edge case values for abstol and span with ITP (#749)

Co-authored-by: Oscar Smith <[email protected]>

Author: dcourteville

lib/BracketingNonlinearSolve/src/itp.jl

@@ -96,8 +96,13 @@ function SciMLBase.__solve(
     span = right - left
     k1 = alg.scaled_k1 * span^(1 - k2) # k1 > 0
     n0 = alg.n0
-    n_h = exponent(span / (2 * ϵ))
-    ϵ_s = ϵ * exp2(n_h + n0)
+    if span / 2 > ϵ * floatmax(typeof(span))
+        # Workaround for when span / (2 * ϵ) == Inf
+        ϵ_s = span / 2 * exp2(n0)
+    else
+        n_h = exponent(span / (2 * ϵ))
+        ϵ_s = ϵ * exp2(n_h) * exp2(n0)
+    end
     T0 = zero(fl)
 
     i = 1

lib/BracketingNonlinearSolve/test/rootfind_tests.jl

@@ -74,6 +74,16 @@ end
             @test result_tol > ϵ
         end
     end
+
+    # Some solvers support abstol=0.0 and converge to floating point precision
+    @testset for alg in (Bisection(), ITP(), Brent())
+        sol = solve(prob, alg; abstol = 0.0)
+        # Test that solution is to floating point precision
+        @test sol.retcode == ReturnCode.FloatingPointLimit
+        @test quadratic_f(sol.left, 2.0) < 0
+        @test quadratic_f(sol.right, 2.0) > 0
+        @test nextfloat(sol.left) == sol.right
+    end
 end
 
 @testitem "Flipped Signs and Reversed Tspan" setup=[RootfindingTestSnippet] tags=[:core] begin