Improve solve_cubic accuracy when the cubic coefficient is negligible

[signalwerk]

Shape::winding/contains misclassify points when a cubic's leading coefficient is floating-point noise

First, a disclosure: I'm a consumer of kurbo (building boolean path operations on top of it) and don't claim to understand its internals. I noticed wrong classifications while testing boolean operations; this report, the analysis, and the suggested fix were crafted with AI assistance (Claude Fable 5). Please treat the fix and the accompanying PR as a suggestion only – feel free to disregard them and solve this however you see fit.

Reproduction

Fresh cargo new with kurbo = "0.13":

use kurbo::{BezPath, Point, Shape};

fn main() {
    // A kappa circle, center (256, 256), r = 204.8, with its on-curve
    // points on the diagonals – an ordinary circle, rotated 45°.
    let d = "M 400.8155 400.8155 \
             C 320.8361 480.7948 191.1639 480.7948 111.1845 400.8155 \
             C 31.2052 320.8361 31.2052 191.1639 111.1845 111.1845 \
             C 191.1639 31.2052 320.8361 31.2052 400.8155 111.1845 \
             C 480.7948 191.1639 480.7948 320.8361 400.8155 400.8155 Z";
    let path = BezPath::from_svg(d).unwrap();

    // p is 204.2676 from the center – clearly inside the circle.
    let p = Point::new(332.85, 445.26);
    assert_ne!(path.winding(p), 0, "winding says p is outside"); // fails
}

Confirmed in 0.12.0, 0.13.1 and current main; silently wrong (no panic). It is not one point: the whole band x ∈ [329, 334] at y = 445.26 is misclassified – the true boundary at that height is at x = 334.394 (by bisection), 6 units away from where kurbo places it. The shape is borrowed from paper.js's boolean-operations test suite¹ (full credit to the paper.js project!).

Cause (as far as we understand it)

The rotation makes each segment's y(t) nearly antisymmetric about t = 0.5, so its cubic coefficient mathematically cancels – but rounding leaves noise. For the deciding y-monotone piece, winding calls

solve_cubic(-44.444500000000005, 119.96895000000012,
            -59.984475000000316, -5.684341886080802e-14)
→ [1.4686287109375038, 0.5313712890624962, -1055258043976628.0]

solve_cubic only falls back to solve_quadratic when 1/c3 overflows; here 1/c3 ≈ -1.8e13 is finite, so the general branch scales the other coefficients by it and amplifies their rounding error. The returned t = 0.5313… (the root in [0, 1], which winding trusts) leaves a polynomial residual of 2.37 – it is not a root. Treating the equation as the quadratic it effectively is gives t = 0.4910141…, matching bisection.

Suggested fix

Check c3 relative to the other coefficients before scaling – the counterpart of what solve_quadratic already documents for its nearly-linear case:

if c3.abs() <= 1e-12 * c0.abs().max(c1.abs()).max(c2.abs()) {
    return solve_quadratic(c0, c1, c2).iter().copied().collect();
}

A PR with this change and a regression test is ready. Attached to this issue is solve-cubic-report.html, a self-contained HTML report that shows the problem visually (before/after root tables, winding scans, SVG close-ups) and embeds the program that generates it, together with the steps to reproduce it.

Possibly related: #411 (intersect_line also goes through solve_cubic), #531, #277.

Full report in HTML

solve-cubic-report.html

Footnotes

1. paper.js test and the issue in paper.js ↩