KineticForces - BUGFIX - Decide cubic and quadratic root counts by conditioned tests

src/KineticForces/BounceAveraging.jl

@@ -17,6 +17,9 @@ const MAX_SPLINE_CELLS = 100_000
 const EXTREMUM_MERGE_TOL = 1e-12
 # Cells a hinted cell search may step before falling back to bisection.
 const HINT_STEP_BUDGET = 8
+# Slack on the root-count tests below. Generous by design: a wrongly admitted repeated
+# root is a near-touch at a stationary point, which the monotone-window filter rejects.
+const ROOT_COUNT_TOL = 64 * eps(Float64)
 
 # ============================================================================
 # BounceData struct
@@ -479,21 +482,11 @@ function _surface_b_field(B_vpar)
         push!(bknot, d)                      # S(0) = B at the cell's left knot
         push!(knot, cell.xR)
 
-        # S'(u) = c + 2b·u + 3a·u², u ∈ [0, h)
-        qa, qb, qc = 3a, 2b, c
-        if abs(qa) <= eps(Float64) * max(abs(qb), abs(qc), 1.0)
-            if qb != 0
-                u = -qc / qb
-                (0.0 <= u < h) && push!(theta, cell.xL + u)
-            end
-        else
-            disc = qb^2 - 4 * qa * qc
-            if disc >= 0
-                sq = sqrt(disc)
-                for u in ((-qb - sq) / (2qa), (-qb + sq) / (2qa))
-                    (0.0 <= u < h) && push!(theta, cell.xL + u)
-                end
-            end
+        # S'(u) = c + 2b·u + 3a·u² on [0, h); a double stationary point returns twice and dedupes below.
+        n, u1, u2, _ = _quadratic_real_roots(3a, 2b, c)
+        for (k, u) in ((1, u1), (2, u2))
+            k <= n || break
+            (0.0 <= u < h) && push!(theta, cell.xL + u)
         end
 
         cell.xR >= 1.0 && break
@@ -530,51 +523,58 @@ end
 end
 
 """
-Real roots of `a·u³ + b·u² + c·u + d = 0`, returned as `(count, r1, r2, r3)` with
-unused slots `NaN`. Degenerate leading coefficients fall through to the quadratic and
-linear cases; three distinct real roots use the trigonometric form, which stays well
-conditioned where Cardano's radicals cancel.
+Real roots of `a·u³ + b·u² + c·u + d = 0` as `(count, r1, r2, r3)`, unused slots `NaN`.
+The count comes from a conditioned ratio test rather than the sign of the discriminant,
+which cancels to rounding noise when two roots nearly coincide; three real roots use the
+trigonometric form. A degenerate leading coefficient falls through to `_quadratic_real_roots`.
 """
 function _real_cubic_roots(a::Float64, b::Float64, c::Float64, d::Float64)
-    scale = max(abs(b), abs(c), abs(d), 1.0)
-    if abs(a) <= eps(Float64) * scale
-        if abs(b) <= eps(Float64) * max(abs(c), abs(d), 1.0)
-            c == 0 && return (0, NaN, NaN, NaN)
-            return (1, -d / c, NaN, NaN)
-        end
-        disc = c * c - 4 * b * d
-        disc < 0 && return (0, NaN, NaN, NaN)
-        sq = sqrt(disc)
-        # Cancellation-free quadratic roots (Numerical Recipes §5.6).
-        q = c == 0 ? -0.5 * sq : -0.5 * (c + copysign(sq, c))
-        r1 = q / b
-        r2 = q == 0 ? r1 : d / q
-        return (2, r1, r2, NaN)
-    end
+    abs(a) <= eps(Float64) * max(abs(b), abs(c), abs(d), 1.0) && return _quadratic_real_roots(b, c, d)
 
     B, C, D = b / a, c / a, d / a
     shift = B / 3
     p = C - B * B / 3
     q = 2 * B^3 / 27 - B * C / 3 + D
-    disc = (q / 2)^2 + (p / 3)^3
-
-    if p == 0 && q == 0
-        return (1, -shift, NaN, NaN)               # triple root
-    elseif abs(disc) <= 8 * eps(Float64) * max((q / 2)^2, abs(p / 3)^3)
-        # Repeated root. The discriminant cancels to ~0 here, so the branches below
-        # would lose it: disc > 0 by a rounding step reports only the simple root.
-        t2 = -3q / (2p)
-        return (3, 3q / p - shift, t2 - shift, t2 - shift)
-    elseif disc > 0
-        s = sqrt(disc)
-        t = cbrt(-q / 2 + s) + cbrt(-q / 2 - s)
-        return (1, t - shift, NaN, NaN)
-    else
-        # Three real roots: t_k = 2r·cos(φ − 2πk/3), r = √(−p/3), φ = acos(−q/2r³)/3.
+
+    if p < 0
+        # Three real roots iff |q/2| ≤ r³, r = √(−p/3): the ratio test stays conditioned where
+        # the discriminant cancels. A ratio a rounding step past 1 is a repeated root; clamp.
         r = sqrt(-p / 3)
-        φ = acos(clamp(-q / (2 * r^3), -1.0, 1.0)) / 3
-        return (3, 2r * cos(φ) - shift, 2r * cos(φ - 2π / 3) - shift, 2r * cos(φ - 4π / 3) - shift)
+        arg = -q / (2 * r^3)
+        if abs(arg) <= 1 + ROOT_COUNT_TOL
+            φ = acos(clamp(arg, -1.0, 1.0)) / 3
+            return (3, 2r * cos(φ) - shift, 2r * cos(φ - 2π / 3) - shift, 2r * cos(φ - 4π / 3) - shift)
+        end
+    end
+
+    # One real root (p ≥ 0: monotone; p < 0: |arg| > 1). The max guards sqrt against rounding.
+    s = sqrt(max(0.0, (q / 2)^2 + (p / 3)^3))
+    t = cbrt(-q / 2 + s) + cbrt(-q / 2 - s)
+    return (1, t - shift, NaN, NaN)
+end
+
+"""
+Real roots of `a·u² + b·u + c = 0` in the `(count, r1, r2, r3)` shape of `_real_cubic_roots`,
+`r3` always `NaN`. A discriminant within rounding of zero is a double root rather than none;
+a degenerate leading coefficient falls through to the linear case.
+"""
+function _quadratic_real_roots(a::Float64, b::Float64, c::Float64)
+    if abs(a) <= eps(Float64) * max(abs(b), abs(c), 1.0)
+        b == 0 && return (0, NaN, NaN, NaN)
+        return (1, -c / b, NaN, NaN)
+    end
+    disc = b * b - 4 * a * c
+    if abs(disc) <= ROOT_COUNT_TOL * max(b * b, abs(4 * a * c))
+        r = -b / (2a)
+        return (2, r, r, NaN)
     end
+    disc < 0 && return (0, NaN, NaN, NaN)
+    sq = sqrt(disc)
+    # Cancellation-free form (Numerical Recipes §5.6).
+    q = b == 0 ? -0.5 * sq : -0.5 * (b + copysign(sq, b))
+    r1 = q / a
+    r2 = q == 0 ? r1 : c / q
+    return (2, r1, r2, NaN)
 end
 
 """

test/runtests_kinetic.jl

@@ -190,6 +190,80 @@
                     @test abs(evalpoly(u, (d, c, b, a))) < 1e-9 * max(1.0, abs(u)^3)
                 end
             end
+
+            # Three real roots, two nearly coincident: the discriminant cancels to rounding noise
+            # and a branch trusting its sign drops two roots. Built from the roots, so exact.
+            for (x1, x2, x3) in ((0.37, 0.37000001, 0.91), (0.2, 0.8, 0.80000001), (-1.5, 0.5, 0.5 + 1e-7))
+                a = 1.0
+                b = -(x1 + x2 + x3)
+                c = x1 * x2 + x1 * x3 + x2 * x3
+                d = -x1 * x2 * x3
+                n, r1, r2, r3 = KF._real_cubic_roots(a, b, c, d)
+                @test n == 3
+                got = filter(isfinite, [r1, r2, r3][1:min(n, 3)])
+                # A near-double root is only determined to about its own separation, so the
+                # match tolerance is loose; what matters is that nothing was dropped.
+                for x in (x1, x2, x3)
+                    @test any(abs.(got .- x) .< 1e-6)
+                end
+                for u in got
+                    @test abs(evalpoly(u, (d, c, b, a))) < 1e-12
+                end
+            end
+        end
+
+        @testset "near-double quadratic roots" begin
+            # b² − 4ac cancels to rounding noise for a near-double root; a branch trusting its
+            # sign returns no roots at all. Built from the roots, so exact.
+            for (x1, x2) in ((0.5, 0.5 + 1e-9), (0.5, 0.5), (-0.3, -0.3 + 1e-10), (2.0, 2.0 + 1e-8))
+                a = 1.0
+                b = -(x1 + x2)
+                c = x1 * x2
+                n, r1, r2, _ = KF._quadratic_real_roots(a, b, c)
+                @test n == 2
+                for x in (x1, x2)
+                    @test min(abs(r1 - x), abs(r2 - x)) < 1e-6
+                end
+                for u in (r1, r2)
+                    @test abs(evalpoly(u, (c, b, a))) < 1e-12
+                end
+            end
+            # Well-separated roots and the degenerate forms are unaffected.
+            sep = KF._quadratic_real_roots(1.0, -3.0, 2.0)
+            @test sep[1] == 2 && sort([sep[2], sep[3]]) == [1.0, 2.0]
+            @test KF._quadratic_real_roots(0.0, 2.0, -4.0)[1:2] == (1, 2.0)
+            @test KF._quadratic_real_roots(1.0, 0.0, 1.0)[1] == 0
+        end
+
+        @testset "two stationary points in one cell" begin
+            # A coarse grid puts both stationary points of an interval inside one cell, so the
+            # crossing is solved between them, not between knots. Checked to occur before asserting.
+            xs2 = collect(range(0.0, 1.0; length=11))
+            B2 = @. 2.0 + 0.45 * cos(2pi * xs2) + 0.20 * cos(6pi * xs2)
+            B2[end] = B2[1]
+            Bv2 = cubic_interp(xs2, B2; bc=PeriodicBC())
+            bf2 = KF._surface_b_field(Bv2)
+            cells = [KF._cell_index(bf2.knot, t) for t in bf2.theta]
+            shared = findfirst(i -> cells[i] == cells[i+1], 1:length(cells)-1)
+            @test shared !== nothing
+            if shared !== nothing
+                i = shared
+                blo, bhi = minmax(bf2.bval[i], bf2.bval[i+1])
+                buf = Float64[]
+                hints = ones(Int, length(bf2.theta) + 1)
+                for btarget in range(blo, bhi; length=400)[2:end-1]
+                    lmda = bo / btarget
+                    roots = copy(KF._bounce_points_at_lambda!(buf, hints, bf2, lmda, bo))
+                    ref = sort!(KF.Roots.find_zeros(θ -> KF._vpar_from_spline(Bv2, lmda, bo, θ), 0.0, 1.0); rev=true)
+                    @test length(roots) == length(ref)
+                    # The crossing inside the shared cell lies strictly between the two
+                    # stationary points that bound it.
+                    @test any(bf2.theta[i] < r < bf2.theta[i+1] for r in roots)
+                    for r in roots
+                        @test abs(KF._vpar_from_spline(Bv2, lmda, bo, r)) < 1e-14
+                    end
+                end
+            end
         end
 
         @testset "agrees with the adaptive scan across λ" begin