diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index c62ae8518..f6f350e73 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/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 """ diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index d502385f7..4f47db9fd 100644 --- a/test/runtests_kinetic.jl +++ b/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