From e88b5fdb049203b586edb5e129de805ab4f8714f Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Tue, 20 Jan 2026 15:27:43 -0700 Subject: [PATCH 1/5] Faster bezier root finding on [0, 1] --- lib/matplotlib/bezier.py | 114 +++++++++++++++++++++++++++++++++++---- 1 file changed, 104 insertions(+), 10 deletions(-) diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py index 7daca897e74d..f09166960f98 100644 --- a/lib/matplotlib/bezier.py +++ b/lib/matplotlib/bezier.py @@ -9,6 +9,82 @@ import numpy as np from matplotlib import _api +from numpy.polynomial.polynomial import polyval as _polyval + + +def _bisect_root_finder(f, a, b, tol=1e-12, max_iter=64): + """Find root of f in [a, b] using bisection. Assumes sign change exists.""" + fa = f(a) + for _ in range(max_iter): + mid = (a + b) * 0.5 + fm = f(mid) + if abs(fm) < tol or (b - a) < tol: + return mid + if fa * fm < 0: + b = mid + else: + a, fa = mid, fm + return (a + b) * 0.5 + + +def _real_roots_in_01(coeffs): + """ + Find real roots of polynomial in [0, 1] using sampling and bisection. + coeffs in ascending order: c0 + c1*x + c2*x**2 + ... + """ + deg = len(coeffs) - 1 + n_samples = max(8, deg * 2) + ts = np.linspace(0, 1, n_samples) + vals = _polyval(ts, coeffs) + + signs = np.sign(vals) + sign_changes = np.where((signs[:-1] != signs[1:]) & (signs[:-1] != 0))[0] + + roots = [] + + def f(t): + return _polyval(t, coeffs) + + max_iter = 53 # float64 fractional precision for [0, 1] interval + for i in sign_changes: + roots.append(_bisect_root_finder(f, ts[i], ts[i + 1], max_iter=max_iter)) + + # Check endpoints + if abs(vals[0]) < 1e-12: + roots.insert(0, 0.0) + if abs(vals[-1]) < 1e-12 and (not roots or abs(roots[-1] - 1.0) > 1e-10): + roots.append(1.0) + + return np.asarray(roots) + + +def _quadratic_roots_in_01(c0, c1, c2): + """Real roots of c0 + c1*x + c2*x**2 in [0, 1].""" + if abs(c2) < 1e-12: # Linear + if abs(c1) < 1e-12: + return np.array([]) + root = -c0 / c1 + return np.array([root]) if 0 <= root <= 1 else np.array([]) + + disc = c1 * c1 - 4 * c2 * c0 + if disc < 0: + return np.array([]) + + sqrt_disc = np.sqrt(disc) + # Numerically stable quadratic formula + if c1 >= 0: + q = -0.5 * (c1 + sqrt_disc) + else: + q = -0.5 * (c1 - sqrt_disc) + + roots = [] + if abs(q) > 1e-12: + roots.append(c0 / q) + if abs(c2) > 1e-12: + roots.append(q / c2) + + roots = np.asarray(roots) + return roots[(roots >= 0) & (roots <= 1)] @lru_cache(maxsize=16) @@ -309,17 +385,35 @@ def axis_aligned_extrema(self): if n <= 1: return np.array([]), np.array([]) Cj = self.polynomial_coefficients - dCj = np.arange(1, n+1)[:, None] * Cj[1:] - dims = [] - roots = [] + dCj = np.arange(1, n + 1)[:, None] * Cj[1:] + + all_dims = [] + all_roots = [] + for i, pi in enumerate(dCj.T): - r = np.roots(pi[::-1]) - roots.append(r) - dims.append(np.full_like(r, i)) - roots = np.concatenate(roots) - dims = np.concatenate(dims) - in_range = np.isreal(roots) & (roots >= 0) & (roots <= 1) - return dims[in_range], np.real(roots)[in_range] + # Trim trailing near-zeros to get actual degree + deg = len(pi) - 1 + while deg > 0 and abs(pi[deg]) < 1e-12: + deg -= 1 + + if deg == 0: + continue + elif deg == 1: + root = -pi[0] / pi[1] + r = np.array([root]) if 0 <= root <= 1 else np.array([]) + elif deg == 2: + r = _quadratic_roots_in_01(pi[0], pi[1], pi[2]) + else: + r = _real_roots_in_01(pi[:deg + 1]) + + if len(r) > 0: + all_roots.append(r) + all_dims.append(np.full(len(r), i)) + + if not all_roots: + return np.array([]), np.array([]) + + return np.concatenate(all_dims), np.concatenate(all_roots) def split_bezier_intersecting_with_closedpath( From de8fd7286ec3f3d85ed5190236bbaab376724ca8 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Tue, 3 Feb 2026 08:08:46 -0700 Subject: [PATCH 2/5] Break out bezier root finding and test it --- lib/matplotlib/bezier.py | 59 +++++++++++++++++++++-------- lib/matplotlib/tests/test_bezier.py | 52 ++++++++++++++++++++++++- 2 files changed, 94 insertions(+), 17 deletions(-) diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py index f09166960f98..3dc30fe5c6a4 100644 --- a/lib/matplotlib/bezier.py +++ b/lib/matplotlib/bezier.py @@ -27,7 +27,7 @@ def _bisect_root_finder(f, a, b, tol=1e-12, max_iter=64): return (a + b) * 0.5 -def _real_roots_in_01(coeffs): +def _bisected_roots_in_01(coeffs): """ Find real roots of polynomial in [0, 1] using sampling and bisection. coeffs in ascending order: c0 + c1*x + c2*x**2 + ... @@ -87,6 +87,47 @@ def _quadratic_roots_in_01(c0, c1, c2): return roots[(roots >= 0) & (roots <= 1)] +def _real_roots_in_01(coeffs): + """ + Find real roots of a polynomial in the interval [0, 1]. + + This is optimized for finding roots only in [0, 1], which is faster than + computing all roots with `numpy.roots` and filtering. For polynomials of + degree <= 2, closed-form solutions are used. For higher degrees, sampling + and bisection are used. + + Parameters + ---------- + coeffs : array-like + Polynomial coefficients in ascending order: + ``c[0] + c[1]*x + c[2]*x**2 + ...`` + Note this is the opposite convention from `numpy.roots`. + + Returns + ------- + roots : ndarray + Sorted array of real roots in [0, 1]. + """ + coeffs = np.asarray(coeffs, dtype=float) + + # Trim trailing near-zeros to get actual degree + deg = len(coeffs) - 1 + while deg > 0 and abs(coeffs[deg]) < 1e-12: + deg -= 1 + + if deg <= 0: + return np.array([]) + elif deg == 1: + root = -coeffs[0] / coeffs[1] + return np.array([root]) if 0 <= root <= 1 else np.array([]) + elif deg == 2: + roots = _quadratic_roots_in_01(coeffs[0], coeffs[1], coeffs[2]) + else: + roots = _bisected_roots_in_01(coeffs[:deg + 1]) + + return np.sort(roots) + + @lru_cache(maxsize=16) def _get_coeff_matrix(n): """ @@ -391,21 +432,7 @@ def axis_aligned_extrema(self): all_roots = [] for i, pi in enumerate(dCj.T): - # Trim trailing near-zeros to get actual degree - deg = len(pi) - 1 - while deg > 0 and abs(pi[deg]) < 1e-12: - deg -= 1 - - if deg == 0: - continue - elif deg == 1: - root = -pi[0] / pi[1] - r = np.array([root]) if 0 <= root <= 1 else np.array([]) - elif deg == 2: - r = _quadratic_roots_in_01(pi[0], pi[1], pi[2]) - else: - r = _real_roots_in_01(pi[:deg + 1]) - + r = _real_roots_in_01(pi) if len(r) > 0: all_roots.append(r) all_dims.append(np.full(len(r), i)) diff --git a/lib/matplotlib/tests/test_bezier.py b/lib/matplotlib/tests/test_bezier.py index 65e2c616e738..42ae468579ca 100644 --- a/lib/matplotlib/tests/test_bezier.py +++ b/lib/matplotlib/tests/test_bezier.py @@ -2,7 +2,57 @@ Tests specific to the bezier module. """ -from matplotlib.bezier import inside_circle, split_bezier_intersecting_with_closedpath +import pytest +import numpy as np +from numpy.testing import assert_allclose + +from matplotlib.bezier import ( + _real_roots_in_01, inside_circle, split_bezier_intersecting_with_closedpath +) + + +def _np_real_roots_in_01(coeffs): + """Reference implementation using np.roots for comparison.""" + coeffs = np.asarray(coeffs) + # np.roots expects descending order (highest power first) + all_roots = np.roots(coeffs[::-1]) + # Filter to real roots in [0, 1] + real_mask = np.abs(all_roots.imag) < 1e-10 + real_roots = all_roots[real_mask].real + in_range = (real_roots >= -1e-10) & (real_roots <= 1 + 1e-10) + return np.sort(np.clip(real_roots[in_range], 0, 1)) + + +@pytest.mark.parametrize("coeffs, expected", [ + ([-0.5, 1], [0.5]), + ([-2, 1], []), # roots: [2.0], not in [0, 1] + ([0.1875, -1, 1], [0.25, 0.75]), + ([1, -2.5, 1], [0.5]), # roots: [0.5, 2.0], only one in [0, 1] + ([1, 0, 1], []), # roots: [+-i], not real + ([-0.08, 0.66, -1.5, 1], [0.2, 0.5, 0.8]), + ([5], []), + ([0, 0, 0], []), + ([0, -0.5, 1], [0.0, 0.5]), + ([0.5, -1.5, 1], [0.5, 1.0]), +]) +def test_real_roots_in_01_known_cases(coeffs, expected): + """Test _real_roots_in_01 against known values and np.roots reference.""" + result = _real_roots_in_01(coeffs) + np_expected = _np_real_roots_in_01(coeffs) + assert_allclose(result, expected, atol=1e-10) + assert_allclose(result, np_expected, atol=1e-10) + + +@pytest.mark.parametrize("degree", range(1, 11)) +def test_real_roots_in_01_random(degree): + """Test random polynomials against np.roots.""" + rng = np.random.default_rng(seed=0) + coeffs = rng.uniform(-10, 10, size=degree + 1) + result = _real_roots_in_01(coeffs) + expected = _np_real_roots_in_01(coeffs) + assert len(result) == len(expected) + if len(result) > 0: + assert_allclose(result, expected, atol=1e-8) def test_split_bezier_with_large_values(): From 52cd7dde7a5eb0cc29addefc091bc3c3a6fab40f Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Tue, 3 Feb 2026 16:54:55 -0700 Subject: [PATCH 3/5] Fix stub --- lib/matplotlib/bezier.pyi | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/matplotlib/bezier.pyi b/lib/matplotlib/bezier.pyi index ad82b873affd..d50328bba8a3 100644 --- a/lib/matplotlib/bezier.pyi +++ b/lib/matplotlib/bezier.pyi @@ -22,6 +22,7 @@ def get_normal_points( cx: float, cy: float, cos_t: float, sin_t: float, length: float ) -> tuple[float, float, float, float]: ... def split_de_casteljau(beta: ArrayLike, t: float) -> tuple[np.ndarray, np.ndarray]: ... +def _real_roots_in_01(coeffs: ArrayLike) -> np.ndarray: ... def find_bezier_t_intersecting_with_closedpath( bezier_point_at_t: Callable[[float], tuple[float, float]], inside_closedpath: Callable[[tuple[float, float]], bool], From d007209f136e1ef552356a7f4707f517d85004f4 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Thu, 5 Mar 2026 20:31:47 -0700 Subject: [PATCH 4/5] Simplify bezier root finding --- lib/matplotlib/bezier.py | 63 +++++++--------------------------------- 1 file changed, 11 insertions(+), 52 deletions(-) diff --git a/lib/matplotlib/bezier.py b/lib/matplotlib/bezier.py index 3dc30fe5c6a4..e2ffe7747e55 100644 --- a/lib/matplotlib/bezier.py +++ b/lib/matplotlib/bezier.py @@ -9,53 +9,6 @@ import numpy as np from matplotlib import _api -from numpy.polynomial.polynomial import polyval as _polyval - - -def _bisect_root_finder(f, a, b, tol=1e-12, max_iter=64): - """Find root of f in [a, b] using bisection. Assumes sign change exists.""" - fa = f(a) - for _ in range(max_iter): - mid = (a + b) * 0.5 - fm = f(mid) - if abs(fm) < tol or (b - a) < tol: - return mid - if fa * fm < 0: - b = mid - else: - a, fa = mid, fm - return (a + b) * 0.5 - - -def _bisected_roots_in_01(coeffs): - """ - Find real roots of polynomial in [0, 1] using sampling and bisection. - coeffs in ascending order: c0 + c1*x + c2*x**2 + ... - """ - deg = len(coeffs) - 1 - n_samples = max(8, deg * 2) - ts = np.linspace(0, 1, n_samples) - vals = _polyval(ts, coeffs) - - signs = np.sign(vals) - sign_changes = np.where((signs[:-1] != signs[1:]) & (signs[:-1] != 0))[0] - - roots = [] - - def f(t): - return _polyval(t, coeffs) - - max_iter = 53 # float64 fractional precision for [0, 1] interval - for i in sign_changes: - roots.append(_bisect_root_finder(f, ts[i], ts[i + 1], max_iter=max_iter)) - - # Check endpoints - if abs(vals[0]) < 1e-12: - roots.insert(0, 0.0) - if abs(vals[-1]) < 1e-12 and (not roots or abs(roots[-1] - 1.0) > 1e-10): - roots.append(1.0) - - return np.asarray(roots) def _quadratic_roots_in_01(c0, c1, c2): @@ -91,10 +44,10 @@ def _real_roots_in_01(coeffs): """ Find real roots of a polynomial in the interval [0, 1]. - This is optimized for finding roots only in [0, 1], which is faster than - computing all roots with `numpy.roots` and filtering. For polynomials of - degree <= 2, closed-form solutions are used. For higher degrees, sampling - and bisection are used. + For polynomials of degree <= 2, closed-form solutions are used. + For higher degrees, `numpy.roots` is used as a fallback. In practice, + matplotlib only ever uses cubic bezier curves and axis_aligned_extrema() + differentiates, so we only ever find roots for degree <= 2. Parameters ---------- @@ -123,7 +76,13 @@ def _real_roots_in_01(coeffs): elif deg == 2: roots = _quadratic_roots_in_01(coeffs[0], coeffs[1], coeffs[2]) else: - roots = _bisected_roots_in_01(coeffs[:deg + 1]) + # np.roots expects descending order (highest power first) + eps = 1e-10 + all_roots = np.roots(coeffs[deg::-1]) + real_mask = np.abs(all_roots.imag) < eps + real_roots = all_roots[real_mask].real + in_range = (real_roots >= -eps) & (real_roots <= 1 + eps) + roots = np.clip(real_roots[in_range], 0, 1) return np.sort(roots) From 9a256484618172acb93f573e4bb1389eefed7b71 Mon Sep 17 00:00:00 2001 From: Scott Shambaugh Date: Fri, 6 Mar 2026 12:27:37 -0700 Subject: [PATCH 5/5] Update bezier tests --- lib/matplotlib/tests/test_bezier.py | 59 ++++++++++------------------- 1 file changed, 20 insertions(+), 39 deletions(-) diff --git a/lib/matplotlib/tests/test_bezier.py b/lib/matplotlib/tests/test_bezier.py index 42ae468579ca..ad5e5acfe49e 100644 --- a/lib/matplotlib/tests/test_bezier.py +++ b/lib/matplotlib/tests/test_bezier.py @@ -11,48 +11,29 @@ ) -def _np_real_roots_in_01(coeffs): - """Reference implementation using np.roots for comparison.""" - coeffs = np.asarray(coeffs) - # np.roots expects descending order (highest power first) - all_roots = np.roots(coeffs[::-1]) - # Filter to real roots in [0, 1] - real_mask = np.abs(all_roots.imag) < 1e-10 - real_roots = all_roots[real_mask].real - in_range = (real_roots >= -1e-10) & (real_roots <= 1 + 1e-10) - return np.sort(np.clip(real_roots[in_range], 0, 1)) - - -@pytest.mark.parametrize("coeffs, expected", [ - ([-0.5, 1], [0.5]), - ([-2, 1], []), # roots: [2.0], not in [0, 1] - ([0.1875, -1, 1], [0.25, 0.75]), - ([1, -2.5, 1], [0.5]), # roots: [0.5, 2.0], only one in [0, 1] - ([1, 0, 1], []), # roots: [+-i], not real - ([-0.08, 0.66, -1.5, 1], [0.2, 0.5, 0.8]), - ([5], []), - ([0, 0, 0], []), - ([0, -0.5, 1], [0.0, 0.5]), - ([0.5, -1.5, 1], [0.5, 1.0]), +@pytest.mark.parametrize("roots, expected_in_01", [ + ([0.5], [0.5]), + ([0.25, 0.75], [0.25, 0.75]), + ([0.2, 0.5, 0.8], [0.2, 0.5, 0.8]), + ([0.1, 0.2, 0.3, 0.4], [0.1, 0.2, 0.3, 0.4]), + ([0.0, 0.5], [0.0, 0.5]), + ([0.5, 1.0], [0.5, 1.0]), + ([2.0], []), # outside [0, 1] + ([0.5, 2.0], [0.5]), # one in, one out + ([-1j, 1j], []), # complex roots + ([0.5, -1j, 1j], [0.5]), # mix of real and complex + ([0.3, 0.3], [0.3, 0.3]), # repeated root ]) -def test_real_roots_in_01_known_cases(coeffs, expected): - """Test _real_roots_in_01 against known values and np.roots reference.""" - result = _real_roots_in_01(coeffs) - np_expected = _np_real_roots_in_01(coeffs) - assert_allclose(result, expected, atol=1e-10) - assert_allclose(result, np_expected, atol=1e-10) +def test_real_roots_in_01(roots, expected_in_01): + roots = np.array(roots) + coeffs = np.poly(roots)[::-1] # np.poly gives descending, we need ascending + result = _real_roots_in_01(coeffs.real) + assert_allclose(result, expected_in_01, atol=1e-10) -@pytest.mark.parametrize("degree", range(1, 11)) -def test_real_roots_in_01_random(degree): - """Test random polynomials against np.roots.""" - rng = np.random.default_rng(seed=0) - coeffs = rng.uniform(-10, 10, size=degree + 1) - result = _real_roots_in_01(coeffs) - expected = _np_real_roots_in_01(coeffs) - assert len(result) == len(expected) - if len(result) > 0: - assert_allclose(result, expected, atol=1e-8) +@pytest.mark.parametrize("coeffs", [[5], [0, 0, 0]]) +def test_real_roots_in_01_no_roots(coeffs): + assert len(_real_roots_in_01(coeffs)) == 0 def test_split_bezier_with_large_values():