Add float (#464)

* version bump

* add float for conversion at spots

Author: jverzani

Project.toml

@@ -1,6 +1,6 @@
 name = "Roots"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "2.2.5"
+version = "2.2.7"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

src/Bracketing/alefeld_potra_shi.jl

@@ -66,7 +66,7 @@ function init_state(::AbstractAlefeldPotraShi, F, x₀, x₁, fx₀, fx₁; c=no
     end
 
     if c === nothing # need c, fc to be defined if one is
-        c = a < zero(a) < b ? _middle(a, b) : secant_step(a, b, fa, fb)
+        c = float(a < zero(a) < b ? _middle(a, b) : secant_step(a, b, fa, fb))
         fc = first(F(c))
     end
 
@@ -111,7 +111,7 @@ function update_state(
     l=NullTracks(),
 ) where {T,S}
     atol, rtol = options.xabstol, options.xreltol
-    μ, λ = oftype(rtol, 0.5), oftype(rtol, 0.7)
+    μ, λ = oftype(float(rtol), 0.5), oftype(float(rtol), 0.7)
     tols = (; λ=λ, atol=atol, rtol=rtol)
 
     a::T, b::T, d::T, ee::T = o.xn0, o.xn1, o.d, o.ee
@@ -209,7 +209,7 @@ struct A2425{K} <: AbstractAlefeldPotraShi end
 function calculateΔ(::A2425{K}, F::Callable_Function, c₀::T, ps) where {K,T}
     a, b, d, ee = ps.a, ps.b, ps.d, ps.ee
     fa, fb, fd, fee = ps.fa, ps.fb, ps.fd, ps.fee
-    tols = (λ=oftype(ps.rtol, 0.7), atol=ps.atol, rtol=ps.rtol)
+    tols = (λ=oftype(float(ps.rtol), 0.7), atol=ps.atol, rtol=ps.rtol)
 
     c = a
     for k in 1:K
@@ -258,7 +258,7 @@ fncalls_per_step(::A57{K}) where {K} = K - 1
 function calculateΔ(::A57{K}, F::Callable_Function, c₀::T, ps) where {K,T}
     a, b, d, ee = ps.a, ps.b, ps.d, ps.ee
     fa, fb, fd, fee = ps.fa, ps.fb, ps.fd, ps.fee
-    tols = (λ=oftype(ps.rtol, 0.7), atol=ps.atol, rtol=ps.rtol)
+    tols = (λ=oftype(float(ps.rtol), 0.7), atol=ps.atol, rtol=ps.rtol)
     c, fc = a, fa
 
     for k in 1:K

src/Bracketing/brent.jl

@@ -48,8 +48,8 @@ function update_state(
     fa, fb, fc = state.fxn0, state.fxn1, state.fc
 
     # next step depends on points; inverse quadratic
-    s::T = inverse_quadratic_step(a, b, c, fa, fb, fc)
-    (isnan(s) || isinf(s)) && (s = secant_step(a, b, fa, fb))
+    s = float(inverse_quadratic_step(a, b, c, fa, fb, fc))
+    (isnan(s) || isinf(s)) && (s = float(secant_step(a, b, fa, fb)))
 
     # guard step
     u, v = (3a + b) / 4, b

src/Derivative/lith.jl

@@ -237,6 +237,7 @@ function init_state(
     ys₀,
 ) where {S,D,R,T}
     xs, ys = init_lith(L, F, x₁, fx₁, x₀, fx₀, ys₀) # [x₀,x₁,…,xₛ₋₁], ...
+
     # skip unit consideration here, as won't fit within storage of ys
     state = LithBoonkkampIJzermanState{S,D + 1,R,T}(
         xs[end],    # xₙ
@@ -703,17 +704,20 @@ function lmm(::LithBoonkkampIJzerman{2,0}, xs, fs)
     x0, x1 = xs
     f0, f1 = fs
 
-    (f0 * x1 - f1 * x0) / (f0 - f1)
+    (f0 * x1 - f1 * x0) / (f0 - f1) |> float
 end
 
 function lmm(::LithBoonkkampIJzerman{3,0}, xs, fs)
+    xs, fs
     x0, x1, x2 = xs
     f0, f1, f2 = fs
 
     (
         f0^2 * f1 * x2 - f0^2 * f2 * x1 - f0 * f1^2 * x2 + f0 * f2^2 * x1 + f1^2 * f2 * x0 -
         f1 * f2^2 * x0
-    ) / (f0^2 * f1 - f0^2 * f2 - f0 * f1^2 + f0 * f2^2 + f1^2 * f2 - f1 * f2^2)
+    ) / (
+        f0^2 * f1 - f0^2 * f2 - f0 * f1^2 + f0 * f2^2 + f1^2 * f2 - f1 * f2^2
+    ) |> float
 end
 
 function lmm(::LithBoonkkampIJzerman{4,0}, xs, fs)
@@ -745,7 +749,7 @@ function lmm(::LithBoonkkampIJzerman{4,0}, xs, fs)
         f0 * f2^3 * f3^2 - f0 * f2^2 * f3^3 - f1^3 * f2^2 * f3 +
         f1^3 * f2 * f3^2 +
         f1^2 * f2^3 * f3 - f1^2 * f2 * f3^3 - f1 * f2^3 * f3^2 + f1 * f2^2 * f3^3
-    )
+    ) |> float
 end
 
 function lmm(::LithBoonkkampIJzerman{5,0}, xs, fs)
@@ -901,7 +905,7 @@ function lmm(::LithBoonkkampIJzerman{5,0}, xs, fs)
         f1 * f2^4 * f3^2 * f4^3 +
         f1 * f2^3 * f3^4 * f4^2 - f1 * f2^3 * f3^2 * f4^4 - f1 * f2^2 * f3^4 * f4^3 +
         f1 * f2^2 * f3^3 * f4^4
-    )
+    ) |> float
 end
 
 function lmm(::LithBoonkkampIJzerman{6,0}, xs, fs)

src/convergence.jl

@@ -235,6 +235,13 @@ function is_small_Δx(
     return δ ≤ Δₓ
 end
 
+
+isnan_x(M::AbstractBracketingMethod, state) = isnan(state.xn1) || isnan(state.xn0)
+isnan_x(M::AbstractNonBracketingMethod, state) = isnan(state.xn1)
+
+isinf_x(M::AbstractBracketingMethod, state) = isinf(state.xn1) || isinf(state.xn0)
+isinf_x(M::AbstractNonBracketingMethod, state) = isinf(state.xn1)
+
 isnan_f(M::AbstractBracketingMethod, state) = isnan(state.fxn1) || isnan(state.fxn0)
 isnan_f(M::AbstractNonBracketingMethod, state) = isnan(state.fxn1)
 
@@ -271,6 +278,8 @@ In `decide_convergence`, stopped values (and `:x_converged` when `strict=false`)
 function assess_convergence(M::Any, state::AbstractUnivariateZeroState, options)
     # return convergence_flag, boolean
     is_exact_zero_f(M, state, options) && return (:exact_zero, true)
+    isnan_x(M, state) && return (:nan, true)
+    isinf_x(M, state) && return (:inf, true)
     isnan_f(M, state) && return (:nan, true)
     isinf_f(M, state) && return (:inf, true)
     is_approx_zero_f(M, state, options) && return (:f_converged, true)
@@ -330,7 +339,23 @@ function decide_convergence(
     val ∈ (:f_converged, :exact_zero, :converged) && return xn1
 
     ## XXX this could be problematic
-    val == :nan && return xn1
+    if val == :nan
+        # return if Δx small
+        Δₓ = abs(xn1 - xn0)
+        δₐ, δᵣ = options.xabstol, options.xreltol
+        u = min(abs(xn0), abs(xn1))
+        δₓ = max(δₐ, 2 * abs(u) * δᵣ) # needs non-zero δₐ to stop near 0
+        Δₓ ≤ δₓ && return xn1
+
+        # or if abs(fxn0) small
+        ϵₐ, ϵᵣ = options.abstol, options.reltol
+        Δ = max(_unitless(ϵₐ), _unitless(xn0) * ϵᵣ)
+        abs(state.fxn0) ≤ Δ * oneunit(state.fxn0) && return xn0
+
+        # else
+        return nan(T) * xn1
+#        return xn1
+    end
     val == :inf_nan && return xn1
 
     ## stopping is a heuristic, x_converged can mask issues

src/find_zeros.jl

@@ -68,7 +68,8 @@ Base.show(io::IO, alpha::Interval) = print(io, "($(alpha.a), $(alpha.b))")
 
 # check if f(a) is non zero using tolerances max(atol, eps()), rtol
 function _non_zero(fa, a::T, atol, rtol) where {T}
-    abs(fa) >= max(atol, abs(a) * rtol * oneunit(fa) / oneunit(a), oneunit(fa) * eps(T))
+    a, r = atol, abs(a) * rtol * oneunit(fa) / oneunit(a), oneunit(fa) * eps(T)
+    return abs(fa) >= max(promote(a,r)...)
 end
 
 # After splitting by zeros we have intervals (zm, zn) this is used to shrink

src/utils.jl

@@ -111,12 +111,12 @@ function quad_vertex(c, fc, b, fb, a, fa)
 end
 
 ## inverse quadratic
-function inverse_quadratic_step(a::T, b, c, fa, fb, fc) where {T}
+function inverse_quadratic_step(a::T, b, c, fa::S, fb, fc) where {T,S}
     s = zero(T)
     s += a * fb * fc / (fa - fb) / (fa - fc) # quad step
     s += b * fa * fc / (fb - fa) / (fb - fc)
     s += c * fa * fb / (fc - fa) / (fc - fb)
-    s
+    float(s)
 end
 
 ## Different functions for approximating f'(xn)

test/test_find_zero.jl

@@ -36,10 +36,12 @@ struct Order3_Test <: Roots.AbstractSecantMethod end
         @test find_zero(sin, 3.0, m) ≈ pi
         @test find_zero(sin, big(3), m) ≈ pi
         @test find_zero(sin, big(3.0), m) ≈ pi
+        @test find_zero(sin, π, m) ≈ pi
         @test find_zero(x -> x^2 - 2.0f0, 2.0f0, m) ≈ sqrt(2) # issue 421
         @test isnan(solve(ZeroProblem(x -> x^2 + 2, 0.5f0)))
     end
 
+
     ## defaults for method argument
     @test find_zero(sin, 3.0) ≈ pi    # order0()
     @test @inferred(find_zero(sin, (3, 4))) ≈ π   # Bisection()