From a601d0e9dc52399caff805f54d67da7275e7900b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 1 Oct 2026 14:55:32 -0400 Subject: [PATCH 1/2] KineticForces - BUGFIX - Decide cubic and quadratic root counts by conditioned tests _real_cubic_roots classified the root count by the sign of the depressed cubic's discriminant (q/2)^2 + (p/3)^3. When two roots nearly coincide the two terms cancel to the rounding floor and the computed sign is noise: for three genuine real roots with two of them 1e-8 apart it came out positive, the one-real-root branch ran, and two real roots were silently dropped. The bounce-point pipeline was safety-netted -- the lost root left the monotone window empty, the NaN emptied the lambda's root set, and the find_zeros fallback ran -- so no physics answer was wrong, but the helper is generic. Decide the count with the conditioned ratio instead: for p < 0 there are three real roots iff |q/2| <= r^3 with r = sqrt(-p/3), tested as the ratio against 1, where the cancellation never happens. A ratio a rounding step past 1 is a repeated root, which the clamped trigonometric form returns directly, so the explicit repeated-root branch and its ad-hoc discriminant tolerance go away. The same cancellation sat at two more sites -- the quadratic branch for a degenerate leading coefficient, and the stationary-point search in _surface_b_field, where a dropped near-double stationary point would have lost a root pair with no fallback at all -- so both now share one _quadratic_real_roots that treats a discriminant within rounding of zero as a double root rather than none. The tolerance is deliberately generous. A wrongly admitted repeated root is a near-touch at a stationary point, which cannot sit strictly inside the monotone window the cell solve filters on, so it can only lose the residual ranking there, never win it. Tests cover the near-double cubic and quadratic cases, built from their roots so the polynomials have them exactly, and the previously unexercised path where both stationary points bounding an interval fall inside a single spline cell. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- src/KineticForces/BounceAveraging.jl | 112 +++++++++++++++------------ test/runtests_kinetic.jl | 79 +++++++++++++++++++ 2 files changed, 141 insertions(+), 50 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index c62ae8518..c50d5898c 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -17,6 +17,12 @@ 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 conditioned root-count tests in _quadratic_real_roots and _real_cubic_roots. +# Where a discriminant cancels to nothing, a value a rounding step past the boundary is a +# repeated root, not a lost one. Generous by design: a wrongly admitted repeated root is +# a near-touch at a stationary point, which cannot sit strictly inside the monotone window +# _cell_level_root filters on, so it can only lose the residual ranking, never win it. +const ROOT_COUNT_TOL = 64 * eps(Float64) # ============================================================================ # BounceData struct @@ -479,21 +485,12 @@ 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², u ∈ [0, h). A double stationary point comes back twice + # and is merged by the dedupe 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 @@ -531,50 +528,65 @@ 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. +unused slots `NaN`. A degenerate leading coefficient falls through to +`_quadratic_real_roots`. The root count comes from a conditioned ratio test on the +depressed cubic, not from the sign of its discriminant, which cancels to rounding noise +when two roots nearly coincide and then reports one root where there are three. Three +real roots use the trigonometric form, which stays well conditioned where Cardano's +radicals cancel. """ 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³ with r = √(−p/3). Testing that ratio against 1 + # stays conditioned where the discriminant (q/2)² + (p/3)³ cancels to nothing; a + # ratio a rounding step past 1 is a repeated root, which the clamp returns as + # t_k = 2r·cos(φ − 2πk/3), φ = acos(−q/2r³)/3, with two of the three coincident. 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 makes t³ + pt + q monotone, and p < 0 with |arg| > 1 is the + # single-crossing case. Cardano's radicals do not cancel here; the max guards the + # square root against the rounding step that the tolerance above just excluded. + 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` with `r3` always `NaN`. A discriminant within rounding of zero is a +double root rather than none: `b² − 4ac` cancels there and the sign it lands on is noise. +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 """ diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index d502385f7..3b2b2070e 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -190,6 +190,85 @@ @test abs(evalpoly(u, (d, c, b, a))) < 1e-9 * max(1.0, abs(u)^3) end end + + # Three genuine, distinct real roots with two of them nearly coincident. The + # discriminant (q/2)² + (p/3)³ cancels to the rounding floor and can land on + # either sign; a branch that trusts that sign reports one root and drops two. + # Built from the roots so the polynomial has them exactly. + 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 the rounding floor for a near-double root and can land on + # either sign; a branch trusting that sign returns no roots at all. Built from + # the roots so the polynomial has them exactly. + 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 bounding a monotone interval inside + # one spline cell, so the crossing is solved on the window between them rather + # than between knots. The fixture is checked to actually produce the case before + # anything is asserted about it. + 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 From fc68afe90480aee22ceeec0eeb4ea11f013b83fb Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 1 Oct 2026 21:44:00 -0400 Subject: [PATCH 2/2] KineticForces - MINOR - Shorten the root-solver annotations to the two-line standard Comment-only. Six annotations ran three to five lines; each now states what the code does in at most two, with the one load-bearing justification -- why the root-count tolerance can be generous -- kept beside the constant it governs. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- src/KineticForces/BounceAveraging.jl | 38 ++++++++++------------------ test/runtests_kinetic.jl | 17 +++++-------- 2 files changed, 19 insertions(+), 36 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index c50d5898c..f6f350e73 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -17,11 +17,8 @@ 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 conditioned root-count tests in _quadratic_real_roots and _real_cubic_roots. -# Where a discriminant cancels to nothing, a value a rounding step past the boundary is a -# repeated root, not a lost one. Generous by design: a wrongly admitted repeated root is -# a near-touch at a stationary point, which cannot sit strictly inside the monotone window -# _cell_level_root filters on, so it can only lose the residual ranking, never win it. +# 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) # ============================================================================ @@ -485,8 +482,7 @@ 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). A double stationary point comes back twice - # and is merged by the dedupe below. + # 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 @@ -527,13 +523,10 @@ end end """ -Real roots of `a·u³ + b·u² + c·u + d = 0`, returned as `(count, r1, r2, r3)` with -unused slots `NaN`. A degenerate leading coefficient falls through to -`_quadratic_real_roots`. The root count comes from a conditioned ratio test on the -depressed cubic, not from the sign of its discriminant, which cancels to rounding noise -when two roots nearly coincide and then reports one root where there are three. Three -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) abs(a) <= eps(Float64) * max(abs(b), abs(c), abs(d), 1.0) && return _quadratic_real_roots(b, c, d) @@ -544,10 +537,8 @@ function _real_cubic_roots(a::Float64, b::Float64, c::Float64, d::Float64) q = 2 * B^3 / 27 - B * C / 3 + D if p < 0 - # Three real roots iff |q/2| ≤ r³ with r = √(−p/3). Testing that ratio against 1 - # stays conditioned where the discriminant (q/2)² + (p/3)³ cancels to nothing; a - # ratio a rounding step past 1 is a repeated root, which the clamp returns as - # t_k = 2r·cos(φ − 2πk/3), φ = acos(−q/2r³)/3, with two of the three coincident. + # 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) arg = -q / (2 * r^3) if abs(arg) <= 1 + ROOT_COUNT_TOL @@ -556,19 +547,16 @@ function _real_cubic_roots(a::Float64, b::Float64, c::Float64, d::Float64) end end - # One real root: p ≥ 0 makes t³ + pt + q monotone, and p < 0 with |arg| > 1 is the - # single-crossing case. Cardano's radicals do not cancel here; the max guards the - # square root against the rounding step that the tolerance above just excluded. + # 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` with `r3` always `NaN`. A discriminant within rounding of zero is a -double root rather than none: `b² − 4ac` cancels there and the sign it lands on is noise. -A degenerate leading coefficient falls through to the linear case. +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) diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index 3b2b2070e..4f47db9fd 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -191,10 +191,8 @@ end end - # Three genuine, distinct real roots with two of them nearly coincident. The - # discriminant (q/2)² + (p/3)³ cancels to the rounding floor and can land on - # either sign; a branch that trusts that sign reports one root and drops two. - # Built from the roots so the polynomial has them exactly. + # 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) @@ -215,9 +213,8 @@ end @testset "near-double quadratic roots" begin - # b² − 4ac cancels to the rounding floor for a near-double root and can land on - # either sign; a branch trusting that sign returns no roots at all. Built from - # the roots so the polynomial has them exactly. + # 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) @@ -239,10 +236,8 @@ end @testset "two stationary points in one cell" begin - # A coarse grid puts both stationary points bounding a monotone interval inside - # one spline cell, so the crossing is solved on the window between them rather - # than between knots. The fixture is checked to actually produce the case before - # anything is asserted about it. + # 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]