From 89939734233d8ab1afb535d5b3cb62bfc497d7e2 Mon Sep 17 00:00:00 2001
From: VPRamon <vallespuigramon@gmail.com>
Date: Sun, 7 Jun 2026 13:50:30 +0200
Subject: [PATCH 1/5] wip

---
 src/approx/fit.rs |  47 ++++++++-
 src/core/mod.rs   |   6 ++
 src/core/roots.rs | 237 ++++++++++++++++++++++++++++++++++++++++++++++
 src/lib.rs        |   5 +-
 4 files changed, 293 insertions(+), 2 deletions(-)
 create mode 100644 src/core/roots.rs

diff --git a/src/approx/fit.rs b/src/approx/fit.rs
index 9370152..960d183 100644
--- a/src/approx/fit.rs
+++ b/src/approx/fit.rs
@@ -5,7 +5,10 @@
 //! by Chebyshev recurrence in `j`, the work is `O(N²)` multiplies-and-adds
 //! and only `O(N)` calls to `cos` (used to build the node grid once).
 
-use crate::core::{nodes, ChebyScalar, ChebySeries, ChebySeriesOn, Domain, NodeKind};
+use crate::core::{nodes, ChebyError, ChebyScalar, ChebySeries, ChebySeriesOn, Domain, NodeKind};
+
+#[cfg(feature = "alloc")]
+use crate::core::ChebySeriesDyn;
 
 /// Compute Chebyshev coefficients from values sampled at roots of `T_N`.
 ///
@@ -53,6 +56,48 @@ pub fn fit_from_fn<T: ChebyScalar, const N: usize>(f: impl Fn(f64) -> T) -> Cheb
     ChebySeries::new(fit_coeffs(&values))
 }
 
+/// Fit a dynamic Chebyshev series of `degree` on normalized `[-1, 1]`.
+#[cfg(feature = "alloc")]
+pub fn fit_dyn_from_fn(
+    degree: usize,
+    mut f: impl FnMut(f64) -> f64,
+) -> Result<ChebySeriesDyn<f64>, ChebyError> {
+    let n = degree.saturating_add(1).max(2);
+    let nf = n as f64;
+    let mut xs = alloc::vec![0.0_f64; n];
+    let mut values = alloc::vec![0.0_f64; n];
+    for (k, x_slot) in xs.iter_mut().enumerate() {
+        let x = (core::f64::consts::PI * (2.0 * k as f64 + 1.0) / (2.0 * nf)).cos();
+        *x_slot = x;
+        values[k] = f(x);
+    }
+
+    let mut coeffs = alloc::vec![0.0_f64; n];
+    for (k, &v) in values.iter().enumerate() {
+        let x = xs[k];
+        coeffs[0] += v;
+        if n > 1 {
+            coeffs[1] += v * x;
+            let two_x = 2.0 * x;
+            let mut tkm1 = 1.0_f64;
+            let mut tk = x;
+            for coeff in coeffs.iter_mut().take(n).skip(2) {
+                let tkp1 = two_x * tk - tkm1;
+                *coeff += v * tkp1;
+                tkm1 = tk;
+                tk = tkp1;
+            }
+        }
+    }
+    coeffs[0] /= nf;
+    let scale = 2.0 / nf;
+    for coeff in coeffs.iter_mut().skip(1) {
+        *coeff *= scale;
+    }
+
+    ChebySeriesDyn::new(coeffs)
+}
+
 /// Fit a Chebyshev series on a physical domain.
 #[inline]
 pub fn fit_from_fn_on<T, X, const N: usize>(
diff --git a/src/core/mod.rs b/src/core/mod.rs
index 3223efd..780d0ed 100644
--- a/src/core/mod.rs
+++ b/src/core/mod.rs
@@ -28,6 +28,9 @@ pub mod nodes;
 pub mod scalar;
 pub mod series;
 
+#[cfg(feature = "alloc")]
+pub mod roots;
+
 pub use domain::Domain;
 pub use error::ChebyError;
 pub use eval::{evaluate, evaluate_both};
@@ -37,3 +40,6 @@ pub use series::{ChebySeries, ChebySeriesOn};
 
 #[cfg(feature = "alloc")]
 pub use series::{ChebySeriesDyn, ChebySeriesDynOn};
+
+#[cfg(feature = "alloc")]
+pub use roots::RootOptions;
diff --git a/src/core/roots.rs b/src/core/roots.rs
new file mode 100644
index 0000000..11ebe38
--- /dev/null
+++ b/src/core/roots.rs
@@ -0,0 +1,237 @@
+//! Root finding for dynamic Chebyshev series on `[-1, 1]`.
+
+#[cfg(feature = "alloc")]
+use alloc::vec::Vec;
+
+use super::ChebySeriesDyn;
+
+const UNIT_ROOT_TOL: f64 = 1e-13;
+const POLY_ZERO_TOL: f64 = 1e-12;
+const DEFAULT_DEDUPE_EPS: f64 = 1e-10;
+
+/// Options for [`ChebySeriesDyn::roots`].
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct RootOptions {
+    /// Tolerance for accepting a root on `[-1, 1]`.
+    pub unit_tol: f64,
+    /// Residual magnitude treated as zero when evaluating the polynomial.
+    pub zero_tol: f64,
+    /// Minimum separation between returned roots after deduplication.
+    pub dedupe_eps: f64,
+}
+
+impl Default for RootOptions {
+    fn default() -> Self {
+        Self {
+            unit_tol: UNIT_ROOT_TOL,
+            zero_tol: POLY_ZERO_TOL,
+            dedupe_eps: DEFAULT_DEDUPE_EPS,
+        }
+    }
+}
+
+#[cfg(feature = "alloc")]
+impl ChebySeriesDyn<f64> {
+    /// Sum of the absolute values of the last `tail` coefficients.
+    ///
+    /// When `tail` exceeds the series length, all coefficients are used.
+    pub fn tail_norm(&self, tail: usize) -> f64 {
+        let coeffs = self.coeffs();
+        let tail = tail.max(1).min(coeffs.len());
+        coeffs
+            .iter()
+            .rev()
+            .take(tail)
+            .map(|c| c.abs())
+            .sum()
+    }
+
+    /// Find real roots in `[-1, 1]` using derivative segmentation and Brent
+    /// refinement on each sign-change interval.
+    pub fn roots(&self) -> Vec<f64> {
+        self.roots_with(RootOptions::default())
+    }
+
+    /// Find real roots in `[-1, 1]` with explicit tolerances.
+    pub fn roots_with(&self, opts: RootOptions) -> Vec<f64> {
+        let mut roots = self.roots_recursive(0, opts);
+        sort_dedup_f64(&mut roots, opts.dedupe_eps);
+        roots
+    }
+
+    fn roots_recursive(&self, depth: usize, opts: RootOptions) -> Vec<f64> {
+        let coeffs = self.coeffs();
+        if coeffs.len() <= 1 {
+            return Vec::new();
+        }
+        if coeffs.len() == 2 {
+            let a = coeffs[1];
+            if a.abs() <= opts.zero_tol {
+                return Vec::new();
+            }
+            let x = -coeffs[0] / a;
+            if (-1.0 - opts.unit_tol..=1.0 + opts.unit_tol).contains(&x) {
+                return vec![x.clamp(-1.0, 1.0)];
+            }
+            return Vec::new();
+        }
+
+        let mut points = vec![-1.0, 1.0];
+        if depth < coeffs.len() {
+            points.extend(self.derivative().roots_recursive(depth + 1, opts));
+        }
+        sort_dedup_f64(&mut points, opts.dedupe_eps);
+
+        let mut roots = Vec::new();
+        for &x in &points {
+            if self.evaluate(x).abs() <= opts.zero_tol {
+                roots.push(x);
+            }
+        }
+
+        for pair in points.windows(2) {
+            let a = pair[0];
+            let b = pair[1];
+            if b - a <= opts.unit_tol {
+                continue;
+            }
+            let fa = self.evaluate(a);
+            let fb = self.evaluate(b);
+            if fa.abs() <= opts.zero_tol || fb.abs() <= opts.zero_tol {
+                continue;
+            }
+            if fa.signum() * fb.signum() < 0.0 {
+                if let Some(root) =
+                    brent_on_unit(a, b, fa, fb, |x| self.evaluate(x), opts.unit_tol)
+                {
+                    roots.push(root.clamp(-1.0, 1.0));
+                }
+            }
+        }
+
+        roots
+    }
+}
+
+fn sort_dedup_f64(values: &mut Vec<f64>, eps: f64) {
+    values.retain(|v| v.is_finite());
+    values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
+    values.dedup_by(|a, b| (*a - *b).abs() <= eps);
+}
+
+fn brent_on_unit<F>(lo: f64, hi: f64, f_lo: f64, f_hi: f64, mut f: F, tol: f64) -> Option<f64>
+where
+    F: FnMut(f64) -> f64,
+{
+    if !lo.is_finite() || !hi.is_finite() || hi < lo {
+        return None;
+    }
+    if f_lo.abs() <= POLY_ZERO_TOL {
+        return Some(lo);
+    }
+    if f_hi.abs() <= POLY_ZERO_TOL {
+        return Some(hi);
+    }
+    if f_lo.signum() * f_hi.signum() > 0.0 {
+        return None;
+    }
+
+    let mut a = lo;
+    let mut b = hi;
+    let mut fa = f_lo;
+    let mut fb = f_hi;
+    let mut c = a;
+    let mut fc = fa;
+    let mut d = b - a;
+    let mut e = d;
+
+    for _ in 0..100 {
+        if fb.signum() * fc.signum() > 0.0 {
+            c = a;
+            fc = fa;
+            d = b - a;
+            e = d;
+        }
+        if fc.abs() < fb.abs() {
+            a = b;
+            b = c;
+            c = a;
+            fa = fb;
+            fb = fc;
+            fc = fa;
+        }
+
+        let tol1 = 2.0 * f64::EPSILON * b.abs() + tol * 0.5;
+        let xm = 0.5 * (c - b);
+        if xm.abs() <= tol1 || fb.abs() <= POLY_ZERO_TOL {
+            return Some(b);
+        }
+
+        if e.abs() >= tol1 && fa.abs() > fb.abs() {
+            let s = fb / fa;
+            let (mut p, mut q) = if (a - c).abs() <= f64::EPSILON {
+                (2.0 * xm * s, 1.0 - s)
+            } else {
+                let q = fa / fc;
+                let r = fb / fc;
+                (
+                    s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0)),
+                    (q - 1.0) * (r - 1.0) * (s - 1.0),
+                )
+            };
+            if p > 0.0 {
+                q = -q;
+            }
+            p = p.abs();
+
+            let min1 = 3.0 * xm * q - (tol1 * q).abs();
+            let min2 = (e * q).abs();
+            if 2.0 * p < min1.min(min2) {
+                e = d;
+                d = p / q;
+            } else {
+                d = xm;
+                e = d;
+            }
+        } else {
+            d = xm;
+            e = d;
+        }
+
+        a = b;
+        fa = fb;
+        if d.abs() > tol1 {
+            b += d;
+        } else {
+            b += tol1.copysign(xm);
+        }
+        fb = f(b);
+        if !fb.is_finite() {
+            return None;
+        }
+    }
+
+    Some(b)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::approx::fit::fit_dyn_from_fn;
+
+    #[test]
+    fn quadratic_roots_on_unit_interval() {
+        let poly = fit_dyn_from_fn(12, |x| x * x - 0.25).unwrap();
+        let roots = poly.roots();
+        assert_eq!(roots.len(), 2, "{roots:?}");
+        assert!((roots[0] + 0.5).abs() < 1e-10, "{roots:?}");
+        assert!((roots[1] - 0.5).abs() < 1e-10, "{roots:?}");
+    }
+
+    #[test]
+    fn tail_norm_uses_last_coefficients() {
+        let poly = fit_dyn_from_fn(8, |x| 2.0 * x * x - 1.0).unwrap();
+        assert!(poly.tail_norm(4).is_finite());
+        assert!(poly.tail_norm(4) >= 0.0);
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index 151a339..abce270 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -32,8 +32,11 @@ pub use core::{
 #[cfg(feature = "alloc")]
 pub use core::ChebySeriesDyn;
 
+#[cfg(feature = "alloc")]
+pub use core::RootOptions;
+
 #[cfg(feature = "approx")]
-pub use approx::fit::fit_coeffs;
+pub use approx::fit::{fit_coeffs, fit_dyn_from_fn};
 
 #[cfg(feature = "piecewise")]
 pub use piecewise::{ChebySegment, ChebySegmentTable};

From e7acf2d9d528c5d2ddb9d283ebaf7bd7ea8ff72e Mon Sep 17 00:00:00 2001
From: VPRamon <vallespuigramon@gmail.com>
Date: Sun, 7 Jun 2026 14:45:24 +0200
Subject: [PATCH 2/5] wip

---
 CHANGELOG.md       |  29 ++-
 README.md          |  31 +++
 benches/cheby.rs   |  44 +++-
 src/approx/fit.rs  |  81 +++----
 src/core/roots.rs  | 591 ++++++++++++++++++++++++++++++++++++++++-----
 src/core/series.rs |  22 ++
 src/lib.rs         |   8 +-
 tests/roots.rs     | 203 ++++++++++++++++
 8 files changed, 887 insertions(+), 122 deletions(-)
 create mode 100644 tests/roots.rs

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 580c5b0..57b2596 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,20 @@ and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.ht
 
 ## [Unreleased]
 
+### Added
+
+- `fit_dyn_from_fn` and `fit_coeffs_dyn` for dynamic Chebyshev fitting on normalized `[-1, 1]`.
+- `ChebySeriesDyn<f64>::roots`, `roots_with`, and `tail_norm` for root finding on the unit interval via colleague-matrix / Durand–Kerner iteration with Chebyshev-node supplementation for tangency roots.
+- `ChebySeriesDyn<f64>::shifted_constant` and `shift_constant_in_place` for level-set root searches without refitting.
+- `ChebySeriesDynOn<f64, X>::roots` and `roots_with` for domain-mapped roots.
+- `RootOptions` with `validated` and `effective` tolerance sanitization.
+- Comprehensive root-finding tests and benchmarks.
+
+### Changed
+
+- Shared fixed-size and dynamic coefficient fitting through a single `accumulate_coeffs` helper.
+- Root finding converts Chebyshev coefficients to a monic power polynomial, finds all roots simultaneously, and verifies residuals against `RootOptions::zero_tol`.
+
 ## [0.3.0] - 2026-05-18
 
 ### Changed
@@ -28,21 +42,6 @@ and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.ht
   `quadrature`, `spectral`, `serde`, `binary`, `fast-dct`, and `full`.
 - The README, examples, benchmarks, and CI configuration were rewritten to
   match the new crate architecture and release posture.
-
-### Added
-
-- `#![forbid(unsafe_code)]` at the crate root.
-- Core abstractions:
-  `Domain<X>`, `ChebyError`, `ChebySeries<T, N>`, `ChebySeriesDyn<T>`,
-  `ChebySeriesOn<T, X, N>`, `ChebyScalar`, `ChebyTime`,
-  `DifferentiateWith<X>`, and `IntegrateWith<X>`.
-- Chebyshev basis functions in `core::basis`:
-  `t(n, x)` and `u(n, x)`.
-- First-class node families in `core::nodes`:
-  `Roots`, `Extrema`, `Lobatto`, `Gauss`, and `GaussLobatto`.
-- Stable Clenshaw evaluation in `core::eval`:
-  `evaluate`, `evaluate_both`, and the compatibility `evaluate_derivative`.
-- Approximation APIs in `approx`:
   coefficient fitting, function fitting on normalized and typed domains,
   barycentric interpolation, adaptive fitting with `FitReport`, coefficient-tail
   error estimates, and a Remez-style minimax interface.
diff --git a/README.md b/README.md
index 2f0dd20..5b831f7 100644
--- a/README.md
+++ b/README.md
@@ -136,6 +136,37 @@ and segments apply the chain rule and use `qtty` multiplication/division so:
 - velocity over time integrates to position,
 - angular rate over time integrates to angle.
 
+## Root finding
+
+With the `alloc` feature, dynamic series support real root finding on the
+normalized interval `[-1, 1]`:
+
+```rust
+use cheby::{fit_dyn_from_fn, RootOptions};
+
+let series = fit_dyn_from_fn(12, |tau| tau * tau - 0.25).unwrap();
+let roots = series.roots(); // normalized tau in [-1, 1]
+
+let opts = RootOptions {
+    zero_tol: 1e-10,
+    ..RootOptions::default()
+};
+let refined = series.roots_with(opts);
+```
+
+The solver converts the Chebyshev expansion to a monic power polynomial, finds
+all roots with a Durand–Kerner iteration (equivalent to colleague-matrix
+eigenvalues for moderate degrees), then verifies residuals on `[-1, 1]`. A
+Chebyshev-node scan supplements tangency and repeated roots when the primary
+pass misses them. [`RootOptions`] controls bracket width, residual tolerance,
+and deduplication. Roots are numerical approximations intended for moderate
+polynomial degrees.
+
+Map normalized roots to a physical domain with [`ChebySeriesDynOn::roots`],
+or shift the constant term without refitting via
+[`ChebySeriesDyn::shifted_constant`] when searching multiple level sets of the
+same fitted signal.
+
 ## Quadrature
 
 `quadrature` contains Chebyshev-related rules only: Clenshaw-Curtis style
diff --git a/benches/cheby.rs b/benches/cheby.rs
index bb24433..8e2f6fa 100644
--- a/benches/cheby.rs
+++ b/benches/cheby.rs
@@ -7,7 +7,7 @@
 
 use std::hint::black_box;
 
-use criterion::{criterion_group, criterion_main, Criterion};
+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
 use qtty::{Kilometer, Second};
 
 fn evaluate_f64(c: &mut Criterion) {
@@ -173,6 +173,44 @@ fn spectral_diff_matrix(c: &mut Criterion) {
     });
 }
 
+fn fit_dyn(c: &mut Criterion) {
+    let mut group = c.benchmark_group("fit_dyn_from_fn_sin");
+    for degree in [4usize, 8, 12, 16, 24, 32] {
+        group.bench_with_input(BenchmarkId::new("deg", degree), &degree, |b, &degree| {
+            b.iter(|| cheby::fit_dyn_from_fn(degree, |x| black_box(x).sin()).unwrap())
+        });
+    }
+    group.finish();
+}
+
+fn evaluate_dyn(c: &mut Criterion) {
+    let series = cheby::fit_dyn_from_fn(16, f64::sin).unwrap();
+    c.bench_function("evaluate_dyn_n17", |b| {
+        b.iter(|| black_box(series.evaluate(black_box(0.37))))
+    });
+}
+
+fn roots_known_quadratic(c: &mut Criterion) {
+    let series = cheby::fit_dyn_from_fn(12, |x| x * x - 0.25).unwrap();
+    c.bench_function("roots_quadratic_deg12", |b| {
+        b.iter(|| black_box(series.roots()))
+    });
+}
+
+fn roots_by_degree(c: &mut Criterion) {
+    let mut group = c.benchmark_group("roots_product");
+    for degree in [4usize, 8, 12, 16, 24, 32] {
+        let series = cheby::fit_dyn_from_fn(degree, |x| {
+            (1..=degree).fold(x, |acc, k| acc * (x - (k as f64 * 0.07 - 0.5)))
+        })
+        .unwrap();
+        group.bench_with_input(BenchmarkId::new("deg", degree), &series, |b, series| {
+            b.iter(|| black_box(series.roots()))
+        });
+    }
+    group.finish();
+}
+
 criterion_group!(
     benches,
     evaluate_f64,
@@ -181,6 +219,10 @@ criterion_group!(
     evaluate_derivative_quantity,
     fit_fixed,
     fit_quantity,
+    fit_dyn,
+    evaluate_dyn,
+    roots_known_quadratic,
+    roots_by_degree,
     adaptive_fit,
     minimax_remez,
     piecewise_lookup,
diff --git a/src/approx/fit.rs b/src/approx/fit.rs
index 960d183..a2a3703 100644
--- a/src/approx/fit.rs
+++ b/src/approx/fit.rs
@@ -16,22 +16,40 @@ use crate::core::ChebySeriesDyn;
 /// sample `k` in a single pass over `j`, avoiding `N²` `cos()` calls.
 #[inline]
 pub fn fit_coeffs<T: ChebyScalar, const N: usize>(values: &[T; N]) -> [T; N] {
+    let xs = nodes::<N>(NodeKind::Roots);
     let mut coeffs = [T::zero(); N];
-    if N == 0 {
-        return coeffs;
+    accumulate_coeffs(values.as_slice(), &xs, &mut coeffs);
+    coeffs
+}
+
+/// Compute dynamic Chebyshev coefficients from samples at supplied unit nodes.
+#[cfg(feature = "alloc")]
+pub fn fit_coeffs_dyn(values: &[f64], unit_nodes: &[f64]) -> Vec<f64> {
+    let mut coeffs = alloc::vec![0.0_f64; values.len()];
+    accumulate_coeffs(values, unit_nodes, &mut coeffs);
+    coeffs
+}
+
+fn accumulate_coeffs<T: ChebyScalar>(values: &[T], unit_nodes: &[f64], coeffs: &mut [T]) {
+    let n = values.len();
+    if n == 0 || coeffs.len() != n {
+        return;
     }
-    let xs = nodes::<N>(NodeKind::Roots);
+    debug_assert_eq!(n, unit_nodes.len());
+
+    for coeff in coeffs.iter_mut() {
+        *coeff = T::zero();
+    }
+
     for (k, &v) in values.iter().enumerate() {
-        let x = xs[k];
-        // T_0 = 1
+        let x = unit_nodes[k];
         coeffs[0] = coeffs[0] + v;
-        if N > 1 {
-            // T_1 = x
+        if n > 1 {
             coeffs[1] = coeffs[1] + v * x;
             let two_x = 2.0 * x;
             let mut tkm1 = 1.0_f64;
             let mut tk = x;
-            for coeff in coeffs.iter_mut().take(N).skip(2) {
+            for coeff in coeffs.iter_mut().take(n).skip(2) {
                 let tkp1 = two_x * tk - tkm1;
                 *coeff = *coeff + v * tkp1;
                 tkm1 = tk;
@@ -39,13 +57,13 @@ pub fn fit_coeffs<T: ChebyScalar, const N: usize>(values: &[T; N]) -> [T; N] {
             }
         }
     }
-    let nf = N as f64;
+
+    let nf = n as f64;
     coeffs[0] = coeffs[0] / nf;
     let scale = 2.0 / nf;
     for coeff in coeffs.iter_mut().skip(1) {
         *coeff = *coeff * scale;
     }
-    coeffs
 }
 
 /// Fit a Chebyshev series on normalized `[-1, 1]`.
@@ -63,41 +81,20 @@ pub fn fit_dyn_from_fn(
     mut f: impl FnMut(f64) -> f64,
 ) -> Result<ChebySeriesDyn<f64>, ChebyError> {
     let n = degree.saturating_add(1).max(2);
-    let nf = n as f64;
-    let mut xs = alloc::vec![0.0_f64; n];
-    let mut values = alloc::vec![0.0_f64; n];
-    for (k, x_slot) in xs.iter_mut().enumerate() {
-        let x = (core::f64::consts::PI * (2.0 * k as f64 + 1.0) / (2.0 * nf)).cos();
-        *x_slot = x;
-        values[k] = f(x);
-    }
-
-    let mut coeffs = alloc::vec![0.0_f64; n];
-    for (k, &v) in values.iter().enumerate() {
-        let x = xs[k];
-        coeffs[0] += v;
-        if n > 1 {
-            coeffs[1] += v * x;
-            let two_x = 2.0 * x;
-            let mut tkm1 = 1.0_f64;
-            let mut tk = x;
-            for coeff in coeffs.iter_mut().take(n).skip(2) {
-                let tkp1 = two_x * tk - tkm1;
-                *coeff += v * tkp1;
-                tkm1 = tk;
-                tk = tkp1;
-            }
-        }
-    }
-    coeffs[0] /= nf;
-    let scale = 2.0 / nf;
-    for coeff in coeffs.iter_mut().skip(1) {
-        *coeff *= scale;
-    }
-
+    let unit_nodes = unit_root_nodes(n);
+    let values: Vec<f64> = unit_nodes.iter().map(|&x| f(x)).collect();
+    let coeffs = fit_coeffs_dyn(&values, &unit_nodes);
     ChebySeriesDyn::new(coeffs)
 }
 
+#[cfg(feature = "alloc")]
+fn unit_root_nodes(n: usize) -> Vec<f64> {
+    let nf = n as f64;
+    (0..n)
+        .map(|k| (core::f64::consts::PI * (2.0 * k as f64 + 1.0) / (2.0 * nf)).cos())
+        .collect()
+}
+
 /// Fit a Chebyshev series on a physical domain.
 #[inline]
 pub fn fit_from_fn_on<T, X, const N: usize>(
diff --git a/src/core/roots.rs b/src/core/roots.rs
index 11ebe38..2ebab97 100644
--- a/src/core/roots.rs
+++ b/src/core/roots.rs
@@ -1,20 +1,58 @@
-//! Root finding for dynamic Chebyshev series on `[-1, 1]`.
+//! Root finding for dynamic Chebyshev series on the normalized interval `[-1, 1]`.
+//!
+//! # Algorithm
+//!
+//! Real roots of a Chebyshev expansion `p(τ) = Σ aₖ Tₖ(τ)` on `τ ∈ [-1, 1]` are
+//! found via the **colleague-matrix / eigenvalue formulation**:
+//!
+//! 1. Trim trailing near-zero high coefficients and reject non-finite input.
+//! 2. Convert the Chebyshev expansion to a monomial (power-basis) polynomial on
+//!    `[-1, 1]`. This avoids evaluating ill-conditioned high-degree monomials directly
+//!    while still enabling robust all-roots solvers.
+//! 3. Find all roots simultaneously with a Durand–Kerner iteration on the monic
+//!    power polynomial (equivalent to colleague-matrix eigenvalues for moderate
+//!    degrees).
+//! 4. Keep real roots in `[-1, 1]`, deduplicate, and verify each against
+//!    `|p(τ)| ≤ zero_tol`.
+//! 5. Supplement with Chebyshev-node scanning and Brent refinement on any missed
+//!    sign-change intervals or tangency minima (important for repeated roots).
+//!
+//! Linear series (`len == 2`) use a closed-form root. Endpoint and repeated roots
+//! are only resolved to the extent allowed by [`RootOptions::zero_tol`] and
+//! [`RootOptions::dedupe_eps`].
+//!
+//! # Coordinate system
+//!
+//! [`ChebySeriesDyn::roots`] and [`ChebySeriesDyn::roots_with`] return roots in the
+//! **normalized** coordinate `τ ∈ [-1, 1]`. Map to a physical domain with
+//! [`Domain::denormalize`](crate::Domain::denormalize) or use [`ChebySeriesDynOn::roots`].
+//!
+//! # Limitations
+//!
+//! - Intended for **moderate degrees** (roughly up to a few dozen coefficients).
+//!   Cost grows superlinearly with degree because all roots are found at once.
+//! - Roots are **numerical approximations** controlled by [`RootOptions`].
+//! - Multiple or tangent roots are harder; tighten `zero_tol` or reduce degree /
+//!   segment length when residuals are unsatisfactory.
+//! - Non-finite coefficients yield an empty root list.
 
 #[cfg(feature = "alloc")]
 use alloc::vec::Vec;
 
-use super::ChebySeriesDyn;
+use super::{ChebyError, ChebySeriesDyn, ChebySeriesDynOn, ChebyTime};
 
-const UNIT_ROOT_TOL: f64 = 1e-13;
-const POLY_ZERO_TOL: f64 = 1e-12;
+const DEFAULT_UNIT_TOL: f64 = 1e-13;
+const DEFAULT_ZERO_TOL: f64 = 1e-12;
 const DEFAULT_DEDUPE_EPS: f64 = 1e-10;
+const MIN_POSITIVE_TOL: f64 = f64::EPSILON;
+const DURAND_KERNER_MAX_ITER: usize = 80;
 
-/// Options for [`ChebySeriesDyn::roots`].
+/// Options for [`ChebySeriesDyn::roots_with`].
 #[derive(Debug, Clone, Copy, PartialEq)]
 pub struct RootOptions {
-    /// Tolerance for accepting a root on `[-1, 1]`.
+    /// Bracket width and root-polishing scale on `[-1, 1]`.
     pub unit_tol: f64,
-    /// Residual magnitude treated as zero when evaluating the polynomial.
+    /// Residual magnitude `|p(τ)|` treated as a root.
     pub zero_tol: f64,
     /// Minimum separation between returned roots after deduplication.
     pub dedupe_eps: f64,
@@ -23,13 +61,58 @@ pub struct RootOptions {
 impl Default for RootOptions {
     fn default() -> Self {
         Self {
-            unit_tol: UNIT_ROOT_TOL,
-            zero_tol: POLY_ZERO_TOL,
+            unit_tol: DEFAULT_UNIT_TOL,
+            zero_tol: DEFAULT_ZERO_TOL,
             dedupe_eps: DEFAULT_DEDUPE_EPS,
         }
     }
 }
 
+impl RootOptions {
+    /// Validate tolerances, returning an error for non-finite, zero, or negative values.
+    pub fn validated(self) -> Result<Self, ChebyError> {
+        Ok(Self {
+            unit_tol: positive_finite_tol(self.unit_tol)?,
+            zero_tol: positive_finite_tol(self.zero_tol)?,
+            dedupe_eps: positive_finite_tol(self.dedupe_eps)?,
+        })
+    }
+
+    /// Return sanitized tolerances, substituting defaults for invalid fields.
+    #[inline]
+    pub fn effective(self) -> Self {
+        self.validated().unwrap_or_default()
+    }
+}
+
+fn positive_finite_tol(value: f64) -> Result<f64, ChebyError> {
+    if !value.is_finite() || value <= 0.0 {
+        Err(ChebyError::NonFiniteInput)
+    } else {
+        Ok(value.max(MIN_POSITIVE_TOL))
+    }
+}
+
+fn coefficients_are_finite(coeffs: &[f64]) -> bool {
+    coeffs.iter().all(|c| c.is_finite())
+}
+
+fn on_unit_interval(x: f64, unit_tol: f64) -> bool {
+    x.is_finite() && (-1.0 - unit_tol..=1.0 + unit_tol).contains(&x)
+}
+
+fn clamp_unit(x: f64) -> f64 {
+    x.clamp(-1.0, 1.0)
+}
+
+fn trim_trailing_coeffs(coeffs: &[f64], zero_tol: f64) -> &[f64] {
+    let mut end = coeffs.len();
+    while end > 1 && coeffs[end - 1].abs() <= zero_tol {
+        end -= 1;
+    }
+    &coeffs[..end]
+}
+
 #[cfg(feature = "alloc")]
 impl ChebySeriesDyn<f64> {
     /// Sum of the absolute values of the last `tail` coefficients.
@@ -38,78 +121,378 @@ impl ChebySeriesDyn<f64> {
     pub fn tail_norm(&self, tail: usize) -> f64 {
         let coeffs = self.coeffs();
         let tail = tail.max(1).min(coeffs.len());
-        coeffs
-            .iter()
-            .rev()
-            .take(tail)
-            .map(|c| c.abs())
-            .sum()
+        coeffs.iter().rev().take(tail).map(|c| c.abs()).sum()
     }
 
-    /// Find real roots in `[-1, 1]` using derivative segmentation and Brent
-    /// refinement on each sign-change interval.
+    /// Find real roots in the normalized interval `[-1, 1]`.
+    ///
+    /// Returns roots as normalized coordinates `τ`, not physical domain values.
+    /// See [`ChebySeriesDynOn::roots`] for domain-mapped roots.
     pub fn roots(&self) -> Vec<f64> {
         self.roots_with(RootOptions::default())
     }
 
     /// Find real roots in `[-1, 1]` with explicit tolerances.
+    ///
+    /// Invalid [`RootOptions`] fields are replaced with defaults via
+    /// [`RootOptions::effective`]. Non-finite coefficients yield an empty list.
     pub fn roots_with(&self, opts: RootOptions) -> Vec<f64> {
-        let mut roots = self.roots_recursive(0, opts);
+        let opts = opts.effective();
+        if !coefficients_are_finite(self.coeffs()) {
+            return Vec::new();
+        }
+
+        let coeffs = trim_trailing_coeffs(self.coeffs(), opts.zero_tol);
+        if coeffs.is_empty() {
+            return Vec::new();
+        }
+        if coeffs.len() == 1 {
+            if coeffs[0].abs() <= opts.zero_tol && on_unit_interval(0.0, opts.unit_tol) {
+                return vec![0.0];
+            }
+            return Vec::new();
+        }
+        if coeffs.len() == 2 {
+            return self.linear_zero(coeffs, opts);
+        }
+
+        let power = chebyshev_to_power(coeffs);
+        let mut roots = durand_kerner_real_roots(&power, opts);
+        let expected = coeffs.len() - 1;
+        if roots.len() < expected {
+            supplement_roots_on_chebyshev_nodes(self, &mut roots, opts);
+        }
         sort_dedup_f64(&mut roots, opts.dedupe_eps);
+        roots.retain(|&r| self.evaluate(r).abs() <= opts.zero_tol);
+        if roots.is_empty() {
+            supplement_roots_on_chebyshev_nodes(self, &mut roots, opts);
+            sort_dedup_f64(&mut roots, opts.dedupe_eps);
+            roots.retain(|&r| self.evaluate(r).abs() <= opts.zero_tol);
+        }
         roots
     }
 
-    fn roots_recursive(&self, depth: usize, opts: RootOptions) -> Vec<f64> {
-        let coeffs = self.coeffs();
-        if coeffs.len() <= 1 {
+    fn linear_zero(&self, coeffs: &[f64], opts: RootOptions) -> Vec<f64> {
+        debug_assert_eq!(coeffs.len(), 2);
+        let a = coeffs[1];
+        if a.abs() <= opts.zero_tol {
+            if coeffs[0].abs() <= opts.zero_tol && on_unit_interval(0.0, opts.unit_tol) {
+                return vec![0.0];
+            }
             return Vec::new();
         }
-        if coeffs.len() == 2 {
-            let a = coeffs[1];
-            if a.abs() <= opts.zero_tol {
-                return Vec::new();
+        let x = -coeffs[0] / a;
+        if on_unit_interval(x, opts.unit_tol) {
+            let root = clamp_unit(x);
+            if self.evaluate(root).abs() <= opts.zero_tol {
+                return vec![root];
             }
-            let x = -coeffs[0] / a;
-            if (-1.0 - opts.unit_tol..=1.0 + opts.unit_tol).contains(&x) {
-                return vec![x.clamp(-1.0, 1.0)];
-            }
-            return Vec::new();
         }
+        Vec::new()
+    }
+}
 
-        let mut points = vec![-1.0, 1.0];
-        if depth < coeffs.len() {
-            points.extend(self.derivative().roots_recursive(depth + 1, opts));
-        }
-        sort_dedup_f64(&mut points, opts.dedupe_eps);
+#[cfg(feature = "alloc")]
+impl<X: ChebyTime> ChebySeriesDynOn<f64, X> {
+    /// Find real roots mapped into the series domain.
+    pub fn roots(&self) -> Vec<X> {
+        self.roots_with(RootOptions::default())
+    }
+
+    /// Find real roots mapped into the series domain with explicit tolerances.
+    pub fn roots_with(&self, opts: RootOptions) -> Vec<X> {
+        self.series()
+            .roots_with(opts)
+            .into_iter()
+            .map(|tau| self.domain().denormalize(tau))
+            .collect()
+    }
+}
 
-        let mut roots = Vec::new();
-        for &x in &points {
-            if self.evaluate(x).abs() <= opts.zero_tol {
-                roots.push(x);
+fn chebyshev_to_power(coeffs: &[f64]) -> Vec<f64> {
+    let degree = coeffs.len().saturating_sub(1);
+    let mut power = vec![0.0; degree + 1];
+    if degree == 0 {
+        power[0] = coeffs[0];
+        return power;
+    }
+
+    let mut t0 = vec![0.0; degree + 1];
+    let mut t1 = vec![0.0; degree + 1];
+    t0[0] = 1.0;
+    t1[1] = 1.0;
+    accumulate_scaled(&mut power, &t0, coeffs[0]);
+    if degree >= 1 {
+        accumulate_scaled(&mut power, &t1, coeffs[1]);
+    }
+
+    for (_, &ck) in coeffs.iter().enumerate().take(degree + 1).skip(2) {
+        let mut tk = vec![0.0; degree + 1];
+        for i in 0..=degree {
+            if i >= 1 {
+                tk[i] += 2.0 * t1[i - 1];
             }
+            tk[i] -= t0[i];
         }
+        accumulate_scaled(&mut power, &tk, ck);
+        t0 = t1;
+        t1 = tk;
+    }
 
-        for pair in points.windows(2) {
-            let a = pair[0];
-            let b = pair[1];
-            if b - a <= opts.unit_tol {
-                continue;
+    power
+}
+
+fn accumulate_scaled(dst: &mut [f64], src: &[f64], scale: f64) {
+    if scale == 0.0 {
+        return;
+    }
+    for (d, &s) in dst.iter_mut().zip(src.iter()) {
+        *d += scale * s;
+    }
+}
+
+fn durand_kerner_real_roots(power: &[f64], opts: RootOptions) -> Vec<f64> {
+    let degree = power.len().saturating_sub(1);
+    if degree == 0 {
+        return Vec::new();
+    }
+
+    let scale = power[degree];
+    if !scale.is_finite() || scale.abs() <= opts.zero_tol {
+        return Vec::new();
+    }
+
+    let mut monic = power.to_vec();
+    for c in &mut monic {
+        *c /= scale;
+    }
+
+    let mut roots = initial_durand_kerner_guess(degree);
+    for _ in 0..DURAND_KERNER_MAX_ITER {
+        let mut next = roots.clone();
+        for i in 0..degree {
+            let value = eval_power(&monic, roots[i].re, roots[i].im);
+            let mut denom = Complex { re: 1.0, im: 0.0 };
+            for (j, zj) in roots.iter().enumerate().take(degree) {
+                if i != j {
+                    denom = cmul(denom, csub(roots[i], *zj));
+                }
             }
-            let fa = self.evaluate(a);
-            let fb = self.evaluate(b);
-            if fa.abs() <= opts.zero_tol || fb.abs() <= opts.zero_tol {
+            if denom.abs() <= f64::EPSILON {
                 continue;
             }
-            if fa.signum() * fb.signum() < 0.0 {
-                if let Some(root) =
-                    brent_on_unit(a, b, fa, fb, |x| self.evaluate(x), opts.unit_tol)
-                {
-                    roots.push(root.clamp(-1.0, 1.0));
+            next[i] = csub(roots[i], cdiv(value, denom));
+        }
+        roots = next;
+    }
+
+    let mut out = Vec::new();
+    for z in roots {
+        if z.im.abs() > opts.unit_tol {
+            continue;
+        }
+        let x = clamp_unit(z.re);
+        if on_unit_interval(x, opts.unit_tol) {
+            out.push(x);
+        }
+    }
+    out
+}
+
+fn initial_durand_kerner_guess(degree: usize) -> Vec<Complex> {
+    let mut roots = Vec::with_capacity(degree);
+    for k in 0..degree {
+        let angle = core::f64::consts::TAU * (k as f64) / (degree as f64);
+        roots.push(Complex {
+            re: 0.8 * angle.cos(),
+            im: 0.8 * angle.sin(),
+        });
+    }
+    roots
+}
+
+#[derive(Clone, Copy, Debug)]
+struct Complex {
+    re: f64,
+    im: f64,
+}
+
+impl Complex {
+    fn abs(self) -> f64 {
+        (self.re * self.re + self.im * self.im).sqrt()
+    }
+}
+
+fn cadd(a: Complex, b: Complex) -> Complex {
+    Complex {
+        re: a.re + b.re,
+        im: a.im + b.im,
+    }
+}
+
+fn csub(a: Complex, b: Complex) -> Complex {
+    Complex {
+        re: a.re - b.re,
+        im: a.im - b.im,
+    }
+}
+
+fn cmul(a: Complex, b: Complex) -> Complex {
+    Complex {
+        re: a.re * b.re - a.im * b.im,
+        im: a.re * b.im + a.im * b.re,
+    }
+}
+
+fn cdiv(a: Complex, b: Complex) -> Complex {
+    let denom = b.re * b.re + b.im * b.im;
+    Complex {
+        re: (a.re * b.re + a.im * b.im) / denom,
+        im: (a.im * b.re - a.re * b.im) / denom,
+    }
+}
+
+fn eval_power(coeffs: &[f64], re: f64, im: f64) -> Complex {
+    let z = Complex { re, im };
+    let mut acc = Complex {
+        re: coeffs[coeffs.len() - 1],
+        im: 0.0,
+    };
+    for &c in coeffs[..coeffs.len() - 1].iter().rev() {
+        acc = cadd(cmul(acc, z), Complex { re: c, im: 0.0 });
+    }
+    acc
+}
+
+fn supplement_roots_on_chebyshev_nodes(
+    series: &ChebySeriesDyn<f64>,
+    roots: &mut Vec<f64>,
+    opts: RootOptions,
+) {
+    let degree = series.coeffs().len().saturating_sub(1).max(1);
+    let nodes = degree.saturating_mul(8).clamp(64, 1024);
+    let mut prev_x = -1.0;
+    let mut prev_f = series.evaluate(prev_x);
+
+    for k in 0..nodes {
+        let x = if k + 1 == nodes {
+            1.0
+        } else {
+            let theta =
+                core::f64::consts::PI * (2.0 * (nodes - 1 - k) as f64 + 1.0) / (2.0 * nodes as f64);
+            theta.cos()
+        };
+        let fx = series.evaluate(x);
+        if fx.abs() <= opts.zero_tol {
+            if !has_nearby_root(roots, x, opts.dedupe_eps) {
+                roots.push(clamp_unit(x));
+            }
+        } else if fx.abs() <= opts.zero_tol * 10.0 {
+            if let Some(root) = polish_newton_minimum(series, x, opts) {
+                if !has_nearby_root(roots, root, opts.dedupe_eps) {
+                    roots.push(root);
+                }
+            }
+        } else if prev_f.abs() > opts.zero_tol && prev_f.signum() * fx.signum() < 0.0 {
+            if let Some(root) = brent_on_unit(
+                prev_x,
+                x,
+                prev_f,
+                fx,
+                |t| series.evaluate(t),
+                opts.unit_tol,
+                opts.zero_tol,
+            ) {
+                if !has_nearby_root(roots, root, opts.dedupe_eps) {
+                    roots.push(clamp_unit(root));
                 }
             }
         }
+        prev_x = x;
+        prev_f = fx;
+    }
 
-        roots
+    let mut best_x = 0.0;
+    let mut best_abs = f64::INFINITY;
+    for k in 0..nodes {
+        let x = if k + 1 == nodes {
+            1.0
+        } else {
+            let theta =
+                core::f64::consts::PI * (2.0 * (nodes - 1 - k) as f64 + 1.0) / (2.0 * nodes as f64);
+            theta.cos()
+        };
+        let fx = series.evaluate(x).abs();
+        if fx < best_abs {
+            best_abs = fx;
+            best_x = x;
+        }
+    }
+    if best_abs <= opts.zero_tol * 100.0 {
+        let delta = (2.0 / nodes as f64).max(opts.unit_tol * 4.0);
+        let lo = clamp_unit(best_x - delta);
+        let hi = clamp_unit(best_x + delta);
+        if let Some(root) = refine_minimum_magnitude(series, lo, hi, opts) {
+            if !has_nearby_root(roots, root, opts.dedupe_eps) {
+                roots.push(root);
+            }
+        } else if best_abs <= opts.zero_tol && !has_nearby_root(roots, best_x, opts.dedupe_eps) {
+            roots.push(clamp_unit(best_x));
+        }
+    }
+}
+
+fn refine_minimum_magnitude(
+    series: &ChebySeriesDyn<f64>,
+    mut lo: f64,
+    mut hi: f64,
+    opts: RootOptions,
+) -> Option<f64> {
+    for _ in 0..48 {
+        if hi - lo <= opts.unit_tol {
+            break;
+        }
+        let m1 = lo + (hi - lo) / 3.0;
+        let m2 = hi - (hi - lo) / 3.0;
+        if series.evaluate(m1).abs() <= series.evaluate(m2).abs() {
+            hi = m2;
+        } else {
+            lo = m1;
+        }
+    }
+    let x = clamp_unit(0.5 * (lo + hi));
+    if series.evaluate(x).abs() <= opts.zero_tol {
+        Some(x)
+    } else {
+        None
+    }
+}
+
+fn has_nearby_root(roots: &[f64], x: f64, eps: f64) -> bool {
+    roots.iter().any(|r| (*r - x).abs() <= eps)
+}
+
+fn polish_newton_minimum(
+    series: &ChebySeriesDyn<f64>,
+    start: f64,
+    opts: RootOptions,
+) -> Option<f64> {
+    let mut x = clamp_unit(start);
+    for _ in 0..32 {
+        let (value, deriv) = series.evaluate_both(x);
+        if value.abs() <= opts.zero_tol {
+            return Some(clamp_unit(x));
+        }
+        if deriv.abs() <= f64::EPSILON {
+            break;
+        }
+        x = clamp_unit(x - value / deriv);
+    }
+    let residual = series.evaluate(x).abs();
+    if residual <= opts.zero_tol {
+        Some(clamp_unit(x))
+    } else {
+        None
     }
 }
 
@@ -119,18 +502,26 @@ fn sort_dedup_f64(values: &mut Vec<f64>, eps: f64) {
     values.dedup_by(|a, b| (*a - *b).abs() <= eps);
 }
 
-fn brent_on_unit<F>(lo: f64, hi: f64, f_lo: f64, f_hi: f64, mut f: F, tol: f64) -> Option<f64>
+fn brent_on_unit<F>(
+    lo: f64,
+    hi: f64,
+    f_lo: f64,
+    f_hi: f64,
+    mut f: F,
+    unit_tol: f64,
+    zero_tol: f64,
+) -> Option<f64>
 where
     F: FnMut(f64) -> f64,
 {
     if !lo.is_finite() || !hi.is_finite() || hi < lo {
         return None;
     }
-    if f_lo.abs() <= POLY_ZERO_TOL {
-        return Some(lo);
+    if f_lo.abs() <= zero_tol {
+        return verified_root(lo, &mut f, unit_tol, zero_tol);
     }
-    if f_hi.abs() <= POLY_ZERO_TOL {
-        return Some(hi);
+    if f_hi.abs() <= zero_tol {
+        return verified_root(hi, &mut f, unit_tol, zero_tol);
     }
     if f_lo.signum() * f_hi.signum() > 0.0 {
         return None;
@@ -144,8 +535,9 @@ where
     let mut fc = fa;
     let mut d = b - a;
     let mut e = d;
+    const BRENT_MAX_ITER: usize = 100;
 
-    for _ in 0..100 {
+    for _ in 0..BRENT_MAX_ITER {
         if fb.signum() * fc.signum() > 0.0 {
             c = a;
             fc = fa;
@@ -161,10 +553,13 @@ where
             fc = fa;
         }
 
-        let tol1 = 2.0 * f64::EPSILON * b.abs() + tol * 0.5;
+        let tol1 = 2.0 * f64::EPSILON * b.abs() + unit_tol * 0.5;
         let xm = 0.5 * (c - b);
-        if xm.abs() <= tol1 || fb.abs() <= POLY_ZERO_TOL {
-            return Some(b);
+        if xm.abs() <= tol1 {
+            return verified_root(b, &mut f, unit_tol, zero_tol);
+        }
+        if fb.abs() <= zero_tol {
+            return verified_root(b, &mut f, unit_tol, zero_tol);
         }
 
         if e.abs() >= tol1 && fa.abs() > fb.abs() {
@@ -205,13 +600,29 @@ where
         } else {
             b += tol1.copysign(xm);
         }
+        b = clamp_unit(b);
         fb = f(b);
         if !fb.is_finite() {
             return None;
         }
     }
 
-    Some(b)
+    verified_root(b, &mut f, unit_tol, zero_tol)
+}
+
+fn verified_root<F>(x: f64, f: &mut F, unit_tol: f64, zero_tol: f64) -> Option<f64>
+where
+    F: FnMut(f64) -> f64,
+{
+    if !on_unit_interval(x, unit_tol) {
+        return None;
+    }
+    let value = f(x);
+    if value.is_finite() && value.abs() <= zero_tol {
+        Some(clamp_unit(x))
+    } else {
+        None
+    }
 }
 
 #[cfg(test)]
@@ -224,6 +635,7 @@ mod tests {
         let poly = fit_dyn_from_fn(12, |x| x * x - 0.25).unwrap();
         let roots = poly.roots();
         assert_eq!(roots.len(), 2, "{roots:?}");
+        assert_roots_within_tol(&poly, &roots, RootOptions::default());
         assert!((roots[0] + 0.5).abs() < 1e-10, "{roots:?}");
         assert!((roots[1] - 0.5).abs() < 1e-10, "{roots:?}");
     }
@@ -234,4 +646,57 @@ mod tests {
         assert!(poly.tail_norm(4).is_finite());
         assert!(poly.tail_norm(4) >= 0.0);
     }
+
+    #[test]
+    fn shift_constant_moves_evaluation() {
+        let poly = fit_dyn_from_fn(6, |x| x).unwrap();
+        let shifted = poly.shifted_constant(2.0);
+        assert!((shifted.evaluate(0.0) - poly.evaluate(0.0) - 2.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn triple_root_fitted() {
+        let p = fit_dyn_from_fn(16, |x| (x + 0.2).powi(3)).unwrap();
+        let opts = RootOptions {
+            zero_tol: 1e-5,
+            unit_tol: 1e-11,
+            dedupe_eps: 1e-10,
+        };
+        let roots = p.roots_with(opts);
+        assert!(!roots.is_empty(), "{roots:?}");
+        assert!(roots.iter().any(|r| (*r + 0.2).abs() < 5e-3), "{roots:?}");
+    }
+
+    #[test]
+    fn chebyshev_polynomial_has_exact_degree_roots() {
+        for n in [4usize, 8, 12] {
+            let mut coeffs = vec![0.0_f64; n + 1];
+            coeffs[n] = 1.0;
+            let p = ChebySeriesDyn::new(coeffs).unwrap();
+            let opts = RootOptions {
+                zero_tol: 1e-9,
+                unit_tol: 1e-11,
+                dedupe_eps: 1e-10,
+            };
+            let roots = p.roots_with(opts);
+            assert_eq!(roots.len(), n, "T_{n} should have {n} roots, got {roots:?}");
+            assert_roots_within_tol(&p, &roots, opts);
+        }
+    }
+
+    pub(crate) fn assert_roots_within_tol(
+        poly: &ChebySeriesDyn<f64>,
+        roots: &[f64],
+        opts: RootOptions,
+    ) {
+        let opts = opts.effective();
+        for &root in roots {
+            assert!(
+                poly.evaluate(root).abs() <= opts.zero_tol,
+                "residual at {root} = {}",
+                poly.evaluate(root)
+            );
+            assert!(on_unit_interval(root, opts.unit_tol));
+        }
+    }
 }
diff --git a/src/core/series.rs b/src/core/series.rs
index 478057d..7bd1d5b 100644
--- a/src/core/series.rs
+++ b/src/core/series.rs
@@ -286,6 +286,28 @@ impl<T: ChebyScalar> ChebySeriesDyn<T> {
     }
 }
 
+#[cfg(feature = "alloc")]
+impl ChebySeriesDyn<f64> {
+    /// Return a copy with the constant (`T₀`) coefficient shifted by `delta`.
+    ///
+    /// Equivalent to evaluating `p(τ) + delta` without refitting.
+    #[must_use]
+    pub fn shifted_constant(&self, delta: f64) -> Self {
+        let mut coeffs = self.coeffs.clone();
+        if !coeffs.is_empty() && delta.is_finite() {
+            coeffs[0] += delta;
+        }
+        Self { coeffs }
+    }
+
+    /// Shift the constant (`T₀`) coefficient in place by `delta`.
+    pub fn shift_constant_in_place(&mut self, delta: f64) {
+        if !self.coeffs.is_empty() && delta.is_finite() {
+            self.coeffs[0] += delta;
+        }
+    }
+}
+
 /// A dynamic Chebyshev series tied to a physical domain.
 #[cfg(feature = "alloc")]
 #[derive(Debug, Clone, PartialEq)]
diff --git a/src/lib.rs b/src/lib.rs
index abce270..ee7ed5b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -32,11 +32,17 @@ pub use core::{
 #[cfg(feature = "alloc")]
 pub use core::ChebySeriesDyn;
 
+#[cfg(feature = "alloc")]
+pub use core::ChebySeriesDynOn;
+
 #[cfg(feature = "alloc")]
 pub use core::RootOptions;
 
 #[cfg(feature = "approx")]
-pub use approx::fit::{fit_coeffs, fit_dyn_from_fn};
+pub use approx::fit::fit_coeffs;
+
+#[cfg(all(feature = "approx", feature = "alloc"))]
+pub use approx::fit::{fit_coeffs_dyn, fit_dyn_from_fn};
 
 #[cfg(feature = "piecewise")]
 pub use piecewise::{ChebySegment, ChebySegmentTable};
diff --git a/tests/roots.rs b/tests/roots.rs
new file mode 100644
index 0000000..49a3c4a
--- /dev/null
+++ b/tests/roots.rs
@@ -0,0 +1,203 @@
+//! Comprehensive tests for dynamic Chebyshev root finding.
+
+use cheby::approx::fit::fit_dyn_from_fn;
+use cheby::{ChebyError, ChebySeriesDyn, Domain, RootOptions};
+use proptest::prelude::*;
+
+fn opts() -> RootOptions {
+    RootOptions {
+        zero_tol: 1e-9,
+        unit_tol: 1e-11,
+        dedupe_eps: 1e-10,
+    }
+}
+
+fn poly(coeffs: &[f64]) -> ChebySeriesDyn<f64> {
+    ChebySeriesDyn::new(coeffs.to_vec()).unwrap()
+}
+
+fn assert_residuals(poly: &ChebySeriesDyn<f64>, roots: &[f64], zero_tol: f64) {
+    for &root in roots {
+        assert!((-1.0..=1.0).contains(&root), "root out of range: {root}");
+        let residual = poly.evaluate(root).abs();
+        assert!(residual <= zero_tol, "residual {residual} at root {root}");
+    }
+}
+
+#[test]
+fn linear_root() {
+    // p(tau) = tau - 0.25  =>  c0 = -0.25, c1 = 1
+    let p = poly(&[-0.25, 1.0]);
+    let roots = p.roots_with(opts());
+    assert_eq!(roots.len(), 1);
+    assert!((roots[0] - 0.25).abs() < 1e-10);
+    assert_residuals(&p, &roots, opts().zero_tol);
+}
+
+#[test]
+fn root_at_minus_one() {
+    // p(tau) = tau + 1
+    let p = poly(&[1.0, 1.0]);
+    let roots = p.roots_with(opts());
+    assert_eq!(roots.len(), 1);
+    assert!((roots[0] + 1.0).abs() < 1e-10);
+    assert_residuals(&p, &roots, opts().zero_tol);
+}
+
+#[test]
+fn root_at_plus_one() {
+    // p(tau) = 1 - tau
+    let p = poly(&[1.0, -1.0]);
+    let roots = p.roots_with(opts());
+    assert_eq!(roots.len(), 1);
+    assert!((roots[0] - 1.0).abs() < 1e-10);
+    assert_residuals(&p, &roots, opts().zero_tol);
+}
+
+#[test]
+fn no_real_roots_on_interval() {
+    // p(tau) = tau^2 + 1  (via direct monomial-ish fit is poor; use evaluation form)
+    let p = fit_dyn_from_fn(8, |x| x * x + 1.0).unwrap();
+    let roots = p.roots_with(opts());
+    assert!(roots.is_empty(), "{roots:?}");
+}
+
+#[test]
+fn double_root_tangency() {
+    let p = fit_dyn_from_fn(12, |x| (x - 0.3).powi(2)).unwrap();
+    let roots = p.roots_with(RootOptions {
+        zero_tol: 1e-6,
+        ..opts()
+    });
+    assert!(!roots.is_empty());
+    assert!(roots.iter().any(|r| (*r - 0.3).abs() < 1e-3));
+    assert_residuals(&p, &roots, 1e-6);
+}
+
+#[test]
+fn triple_root() {
+    let p = fit_dyn_from_fn(16, |x| (x + 0.2).powi(3)).unwrap();
+    let roots = p.roots_with(RootOptions {
+        zero_tol: 1e-5,
+        ..opts()
+    });
+    assert!(!roots.is_empty());
+    assert!(roots.iter().any(|r| (*r + 0.2).abs() < 5e-3));
+}
+
+#[test]
+fn multiple_separated_roots() {
+    let p = fit_dyn_from_fn(16, |x| x * x * x - x).unwrap();
+    let roots = p.roots_with(opts());
+    assert_eq!(roots.len(), 3, "{roots:?}");
+    assert_residuals(&p, &roots, opts().zero_tol);
+}
+
+#[test]
+fn very_close_roots() {
+    let p = fit_dyn_from_fn(24, |x| (x - 0.40) * (x - 0.41)).unwrap();
+    let roots = p.roots_with(RootOptions {
+        zero_tol: 1e-7,
+        dedupe_eps: 1e-4,
+        ..opts()
+    });
+    assert!(roots.len() >= 2, "{roots:?}");
+    assert_residuals(&p, &roots, 1e-7);
+}
+
+#[test]
+fn chebyshev_polynomial_roots() {
+    for n in [4usize, 8, 12] {
+        let mut coeffs = vec![0.0_f64; n + 1];
+        coeffs[n] = 1.0;
+        let p = ChebySeriesDyn::new(coeffs).unwrap();
+        let options = opts();
+        let roots = p.roots_with(options);
+        assert_eq!(roots.len(), n, "T_{n} should have {n} roots, got {roots:?}");
+        assert_residuals(&p, &roots, options.zero_tol);
+    }
+}
+
+#[test]
+fn non_finite_coefficients_yield_no_roots() {
+    let p = poly(&[1.0, f64::NAN]);
+    assert!(p.roots_with(opts()).is_empty());
+    let p = poly(&[f64::INFINITY, 1.0]);
+    assert!(p.roots_with(opts()).is_empty());
+}
+
+#[test]
+fn invalid_root_options_are_sanitized() {
+    let p = poly(&[-0.25, 1.0]);
+    let bad = RootOptions {
+        unit_tol: f64::NAN,
+        zero_tol: -1.0,
+        dedupe_eps: 0.0,
+    };
+    assert_eq!(bad.validated(), Err(ChebyError::NonFiniteInput));
+    let roots = p.roots_with(bad);
+    assert_eq!(roots.len(), 1);
+}
+
+#[test]
+fn deduplication_merges_nearby_roots() {
+    let p = fit_dyn_from_fn(10, |x| x - 0.5).unwrap();
+    let roots = p.roots_with(RootOptions {
+        dedupe_eps: 1.0,
+        ..opts()
+    });
+    assert_eq!(roots.len(), 1);
+}
+
+#[test]
+fn shifted_constant_finds_offset_roots() {
+    let base = fit_dyn_from_fn(10, |x| x).unwrap();
+    let shifted = base.shifted_constant(-0.5);
+    let roots = shifted.roots_with(opts());
+    assert_eq!(roots.len(), 1);
+    assert!((roots[0] - 0.5).abs() < 1e-8);
+}
+
+#[test]
+fn domain_mapped_roots() {
+    let domain = Domain::new(0.0_f64, 10.0);
+    let series = fit_dyn_from_fn(10, |tau| tau - 0.0).unwrap();
+    let on = cheby::ChebySeriesDynOn::new(domain, series);
+    let roots = on.roots_with(opts());
+    assert_eq!(roots.len(), 1);
+    assert!((roots[0] - 5.0).abs() < 1e-8);
+}
+
+proptest! {
+    #[test]
+    fn fitted_quadratic_roots_satisfy_residual(a in -0.9f64..0.9, b in -0.9f64..0.9) {
+        let p = fit_dyn_from_fn(12, move |x| (x - a) * (x - b)).unwrap();
+        let options = opts();
+        let roots = p.roots_with(options);
+        prop_assert!(!roots.is_empty());
+        for root in roots {
+            prop_assert!((-1.0..=1.0).contains(&root));
+            prop_assert!(p.evaluate(root).abs() <= options.zero_tol);
+        }
+    }
+}
+
+#[test]
+fn fixed_and_dynamic_fit_agree_on_samples() {
+    use cheby::approx::fit::{fit_coeffs, fit_coeffs_dyn, fit_from_fn};
+    let values: [f64; 8] = core::array::from_fn(|k| {
+        let n = 8.0;
+        let x = (std::f64::consts::PI * (2.0 * k as f64 + 1.0) / (2.0 * n)).cos();
+        x.sin()
+    });
+    let xs = cheby::nodes::<8>();
+    let fixed = fit_coeffs(&values);
+    let dynamic = fit_coeffs_dyn(&values, &xs);
+    for (a, b) in fixed.iter().zip(dynamic.iter()) {
+        assert!((a - b).abs() < 1e-14, "fixed={a} dynamic={b}");
+    }
+    let series = fit_from_fn::<f64, 8>(|x| x.sin());
+    for (a, b) in series.coeffs().iter().zip(dynamic.iter()) {
+        assert!((a - b).abs() < 1e-14);
+    }
+}

From bdfcfa6c469b5d642653b7287af08b0b8414c3d3 Mon Sep 17 00:00:00 2001
From: VPRamon <vallespuigramon@gmail.com>
Date: Sun, 7 Jun 2026 14:50:27 +0200
Subject: [PATCH 3/5] release

---
 CHANGELOG.md      |  8 +++---
 Cargo.toml        |  2 +-
 README.md         | 17 +++++++------
 src/approx/fit.rs |  7 ++++--
 src/core/roots.rs | 60 ++++++++++++++++++++++++--------------------
 tests/roots.rs    | 64 +++++++++++++++++++++++++++++++++++++++++++++++
 6 files changed, 118 insertions(+), 40 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 57b2596..c20b015 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,12 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
+## [0.4.0] - 2026-06-07
 
 ### Added
 
 - `fit_dyn_from_fn` and `fit_coeffs_dyn` for dynamic Chebyshev fitting on normalized `[-1, 1]`.
-- `ChebySeriesDyn<f64>::roots`, `roots_with`, and `tail_norm` for root finding on the unit interval via colleague-matrix / Durand–Kerner iteration with Chebyshev-node supplementation for tangency roots.
+- `ChebySeriesDyn<f64>::roots`, `roots_with`, and `tail_norm` for root finding on the unit interval via Durand–Kerner iteration with Chebyshev-node fallback for tangency roots.
 - `ChebySeriesDyn<f64>::shifted_constant` and `shift_constant_in_place` for level-set root searches without refitting.
 - `ChebySeriesDynOn<f64, X>::roots` and `roots_with` for domain-mapped roots.
 - `RootOptions` with `validated` and `effective` tolerance sanitization.
@@ -19,7 +19,9 @@ and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.ht
 ### Changed
 
 - Shared fixed-size and dynamic coefficient fitting through a single `accumulate_coeffs` helper.
-- Root finding converts Chebyshev coefficients to a monic power polynomial, finds all roots simultaneously, and verifies residuals against `RootOptions::zero_tol`.
+- Root finding converts Chebyshev coefficients to power basis, applies Durand–Kerner, verifies residuals against `RootOptions::zero_tol`, and falls back to node scanning when needed.
+- Constant and identically-zero series return an empty root list.
+- `fit_dyn_from_fn(0, …)` now produces a single-coefficient constant fit.
 
 ## [0.3.0] - 2026-05-18
 
diff --git a/Cargo.toml b/Cargo.toml
index 24d5d10..482e00e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "cheby"
-version = "0.3.0"
+version = "0.4.0"
 edition = "2021"
 authors = ["VPRamon <vallespuigramon@gmail.com>"]
 description = "Unit-safe Chebyshev approximation and spectral numerics for Rust."
diff --git a/README.md b/README.md
index 5b831f7..1f165f6 100644
--- a/README.md
+++ b/README.md
@@ -154,13 +154,16 @@ let opts = RootOptions {
 let refined = series.roots_with(opts);
 ```
 
-The solver converts the Chebyshev expansion to a monic power polynomial, finds
-all roots with a Durand–Kerner iteration (equivalent to colleague-matrix
-eigenvalues for moderate degrees), then verifies residuals on `[-1, 1]`. A
-Chebyshev-node scan supplements tangency and repeated roots when the primary
-pass misses them. [`RootOptions`] controls bracket width, residual tolerance,
-and deduplication. Roots are numerical approximations intended for moderate
-polynomial degrees.
+The solver converts the Chebyshev expansion to a monic power polynomial and finds
+roots with a **Durand–Kerner** iteration, then verifies residuals on `[-1, 1]`.
+When the primary pass misses roots (tangency, repeated roots), a fallback
+Chebyshev-node scan uses Brent refinement and minimum search. Constant or
+identically-zero series return an empty root list.
+
+[`RootOptions`] controls bracket width (`unit_tol`), residual tolerance
+(`zero_tol`), and deduplication (`dedupe_eps`). Root finding is intended for
+moderate polynomial degrees and can be less stable than Clenshaw evaluation at
+high degree.
 
 Map normalized roots to a physical domain with [`ChebySeriesDynOn::roots`],
 or shift the constant term without refitting via
diff --git a/src/approx/fit.rs b/src/approx/fit.rs
index a2a3703..64aa687 100644
--- a/src/approx/fit.rs
+++ b/src/approx/fit.rs
@@ -74,13 +74,16 @@ pub fn fit_from_fn<T: ChebyScalar, const N: usize>(f: impl Fn(f64) -> T) -> Cheb
     ChebySeries::new(fit_coeffs(&values))
 }
 
-/// Fit a dynamic Chebyshev series of `degree` on normalized `[-1, 1]`.
+/// Fit a dynamic Chebyshev series of polynomial `degree` on normalized `[-1, 1]`.
+///
+/// The returned series has `degree + 1` coefficients (`T_0` through `T_degree`).
+/// `degree = 0` yields a constant (single-coefficient) fit sampled at `τ = 0`.
 #[cfg(feature = "alloc")]
 pub fn fit_dyn_from_fn(
     degree: usize,
     mut f: impl FnMut(f64) -> f64,
 ) -> Result<ChebySeriesDyn<f64>, ChebyError> {
-    let n = degree.saturating_add(1).max(2);
+    let n = degree.saturating_add(1);
     let unit_nodes = unit_root_nodes(n);
     let values: Vec<f64> = unit_nodes.iter().map(|&x| f(x)).collect();
     let coeffs = fit_coeffs_dyn(&values, &unit_nodes);
diff --git a/src/core/roots.rs b/src/core/roots.rs
index 2ebab97..d1f234d 100644
--- a/src/core/roots.rs
+++ b/src/core/roots.rs
@@ -3,23 +3,27 @@
 //! # Algorithm
 //!
 //! Real roots of a Chebyshev expansion `p(τ) = Σ aₖ Tₖ(τ)` on `τ ∈ [-1, 1]` are
-//! found via the **colleague-matrix / eigenvalue formulation**:
+//! found in two stages:
 //!
-//! 1. Trim trailing near-zero high coefficients and reject non-finite input.
-//! 2. Convert the Chebyshev expansion to a monomial (power-basis) polynomial on
-//!    `[-1, 1]`. This avoids evaluating ill-conditioned high-degree monomials directly
-//!    while still enabling robust all-roots solvers.
-//! 3. Find all roots simultaneously with a Durand–Kerner iteration on the monic
-//!    power polynomial (equivalent to colleague-matrix eigenvalues for moderate
-//!    degrees).
-//! 4. Keep real roots in `[-1, 1]`, deduplicate, and verify each against
-//!    `|p(τ)| ≤ zero_tol`.
-//! 5. Supplement with Chebyshev-node scanning and Brent refinement on any missed
-//!    sign-change intervals or tangency minima (important for repeated roots).
+//! 1. **Primary pass:** trim trailing near-zero high coefficients, reject non-finite
+//!    input, convert the Chebyshev expansion to a monic power-basis polynomial, and
+//!    locate roots with a **Durand–Kerner** iteration. Candidates are kept when they
+//!    are real (`|Im| ≤ unit_tol`), lie in `[-1, 1]`, and satisfy
+//!    `|p(τ)| ≤ zero_tol` after Clenshaw evaluation.
+//! 2. **Fallback pass:** if too few roots are found, scan Chebyshev nodes on `[-1, 1]`,
+//!    refine sign-change brackets with **Brent's method**, and search for tangency
+//!    minima with Newton and golden-section refinement. Fallback thresholds scale
+//!    with [`RootOptions::zero_tol`].
 //!
-//! Linear series (`len == 2`) use a closed-form root. Endpoint and repeated roots
-//! are only resolved to the extent allowed by [`RootOptions::zero_tol`] and
-//! [`RootOptions::dedupe_eps`].
+//! Linear series (`len == 2`) use a closed-form root. Constant series return an
+//! empty list (see below). Repeated or tangent roots are only resolved to the extent
+//! allowed by [`RootOptions::zero_tol`] and [`RootOptions::dedupe_eps`].
+//!
+//! # Constant and degenerate series
+//!
+//! A series with a single coefficient (constant), an identically-zero series, or a
+//! near-zero constant within `zero_tol` has **no isolated roots** on `[-1, 1]` and
+//! returns an empty list. A non-zero constant likewise returns an empty list.
 //!
 //! # Coordinate system
 //!
@@ -31,9 +35,11 @@
 //!
 //! - Intended for **moderate degrees** (roughly up to a few dozen coefficients).
 //!   Cost grows superlinearly with degree because all roots are found at once.
+//! - Root finding converts to power basis internally; this can be **less stable than
+//!   Clenshaw evaluation** at high degree even when residuals pass `zero_tol`.
 //! - Roots are **numerical approximations** controlled by [`RootOptions`].
-//! - Multiple or tangent roots are harder; tighten `zero_tol` or reduce degree /
-//!   segment length when residuals are unsatisfactory.
+//! - Multiple or tangent roots are harder; tighten `zero_tol` or reduce degree when
+//!   residuals are unsatisfactory.
 //! - Non-finite coefficients yield an empty root list.
 
 #[cfg(feature = "alloc")]
@@ -136,6 +142,9 @@ impl ChebySeriesDyn<f64> {
     ///
     /// Invalid [`RootOptions`] fields are replaced with defaults via
     /// [`RootOptions::effective`]. Non-finite coefficients yield an empty list.
+    ///
+    /// Constant, identically-zero, and near-zero constant series return an empty
+    /// list because they have no isolated roots on `[-1, 1]`.
     pub fn roots_with(&self, opts: RootOptions) -> Vec<f64> {
         let opts = opts.effective();
         if !coefficients_are_finite(self.coeffs()) {
@@ -143,13 +152,7 @@ impl ChebySeriesDyn<f64> {
         }
 
         let coeffs = trim_trailing_coeffs(self.coeffs(), opts.zero_tol);
-        if coeffs.is_empty() {
-            return Vec::new();
-        }
-        if coeffs.len() == 1 {
-            if coeffs[0].abs() <= opts.zero_tol && on_unit_interval(0.0, opts.unit_tol) {
-                return vec![0.0];
-            }
+        if coeffs.len() <= 1 {
             return Vec::new();
         }
         if coeffs.len() == 2 {
@@ -176,9 +179,6 @@ impl ChebySeriesDyn<f64> {
         debug_assert_eq!(coeffs.len(), 2);
         let a = coeffs[1];
         if a.abs() <= opts.zero_tol {
-            if coeffs[0].abs() <= opts.zero_tol && on_unit_interval(0.0, opts.unit_tol) {
-                return vec![0.0];
-            }
             return Vec::new();
         }
         let x = -coeffs[0] / a;
@@ -667,6 +667,12 @@ mod tests {
         assert!(roots.iter().any(|r| (*r + 0.2).abs() < 5e-3), "{roots:?}");
     }
 
+    #[test]
+    fn constant_zero_series_has_no_roots() {
+        let p = ChebySeriesDyn::new(vec![0.0]).unwrap();
+        assert!(p.roots_with(RootOptions::default()).is_empty());
+    }
+
     #[test]
     fn chebyshev_polynomial_has_exact_degree_roots() {
         for n in [4usize, 8, 12] {
diff --git a/tests/roots.rs b/tests/roots.rs
index 49a3c4a..d9f3346 100644
--- a/tests/roots.rs
+++ b/tests/roots.rs
@@ -168,6 +168,70 @@ fn domain_mapped_roots() {
     assert!((roots[0] - 5.0).abs() < 1e-8);
 }
 
+#[test]
+fn constant_zero_series_has_no_roots() {
+    assert!(poly(&[0.0]).roots_with(opts()).is_empty());
+    assert!(poly(&[1e-15]).roots_with(opts()).is_empty());
+}
+
+#[test]
+fn constant_nonzero_series_has_no_roots() {
+    assert!(poly(&[1.0]).roots_with(opts()).is_empty());
+    assert!(poly(&[-2.5]).roots_with(opts()).is_empty());
+}
+
+#[test]
+fn degenerate_linear_has_no_roots() {
+    let p = poly(&[1.0, 1e-15]);
+    assert!(p.roots_with(opts()).is_empty());
+}
+
+#[test]
+fn fit_degree_zero_produces_one_coefficient() {
+    let p = fit_dyn_from_fn(0, |_| 3.5).unwrap();
+    assert_eq!(p.coeffs().len(), 1);
+    assert!((p.evaluate(0.0) - 3.5).abs() < 1e-12);
+    assert!(p.roots_with(opts()).is_empty());
+}
+
+#[test]
+fn fit_degree_one_produces_two_coefficients() {
+    let p = fit_dyn_from_fn(1, |x| 2.0 * x + 1.0).unwrap();
+    assert_eq!(p.coeffs().len(), 2);
+    let roots = p.roots_with(opts());
+    assert_eq!(roots.len(), 1);
+    assert!((roots[0] + 0.5).abs() < 1e-8);
+    assert_residuals(&p, &roots, opts().zero_tol);
+}
+
+#[test]
+fn chebyshev_t24_roots() {
+    let n = 24usize;
+    let mut coeffs = vec![0.0_f64; n + 1];
+    coeffs[n] = 1.0;
+    let p = ChebySeriesDyn::new(coeffs).unwrap();
+    let options = opts();
+    let roots = p.roots_with(options);
+    assert_eq!(roots.len(), n, "{roots:?}");
+    assert_residuals(&p, &roots, options.zero_tol);
+}
+
+#[test]
+fn strict_zero_tol_can_exclude_approximate_roots() {
+    let p = fit_dyn_from_fn(12, |x| (x - 0.3).powi(2)).unwrap();
+    let loose = p.roots_with(RootOptions {
+        zero_tol: 1e-6,
+        ..opts()
+    });
+    let strict = p.roots_with(RootOptions {
+        zero_tol: 1e-14,
+        ..opts()
+    });
+    assert!(!loose.is_empty());
+    assert!(strict.is_empty() || strict.len() <= loose.len());
+    assert_residuals(&p, &loose, 1e-6);
+}
+
 proptest! {
     #[test]
     fn fitted_quadratic_roots_satisfy_residual(a in -0.9f64..0.9, b in -0.9f64..0.9) {

From a27c22e3ae430eb4a73a26ec509c81e152736056 Mon Sep 17 00:00:00 2001
From: VPRamon <vallespuigramon@gmail.com>
Date: Sun, 7 Jun 2026 14:51:04 +0200
Subject: [PATCH 4/5] update deps

---
 Cargo.lock | 86 +++++++++++++++++++++++++++---------------------------
 1 file changed, 43 insertions(+), 43 deletions(-)

diff --git a/Cargo.lock b/Cargo.lock
index d35d2ee..d443805 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -49,9 +49,9 @@ dependencies = [
 
 [[package]]
 name = "autocfg"
-version = "1.5.0"
+version = "1.5.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c08606f8c3cbf4ce6ec8e28fb0014a2c086708fe954eaa885384a6165172e7e8"
+checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53"
 
 [[package]]
 name = "bit-set"
@@ -70,15 +70,15 @@ checksum = "5e764a1d40d510daf35e07be9eb06e75770908c27d411ee6c92109c9840eaaf7"
 
 [[package]]
 name = "bitflags"
-version = "2.11.1"
+version = "2.13.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c4512299f36f043ab09a583e57bceb5a5aab7a73db1805848e8fef3c9e8c78b3"
+checksum = "b4388bee8683e3d04af747c73422af53102d2bd24d9eadb6cbc100baef4b43f8"
 
 [[package]]
 name = "bumpalo"
-version = "3.20.2"
+version = "3.20.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5d20789868f4b01b2f2caec9f5c4e0213b41e3e5702a50157d699ae31ced2fcb"
+checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649"
 
 [[package]]
 name = "cast"
@@ -88,9 +88,9 @@ checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5"
 
 [[package]]
 name = "cc"
-version = "1.2.62"
+version = "1.2.63"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a1dce859f0832a7d088c4f1119888ab94ef4b5d6795d1ce05afb7fe159d79f98"
+checksum = "556e016178bb5662a08681bbe0f00f8e17631781a4dfc8c45e466e4b185ec27f"
 dependencies = [
  "find-msvc-tools",
  "shlex",
@@ -104,7 +104,7 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801"
 
 [[package]]
 name = "cheby"
-version = "0.3.0"
+version = "0.4.0"
 dependencies = [
  "approx",
  "criterion",
@@ -234,9 +234,9 @@ checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5"
 
 [[package]]
 name = "either"
-version = "1.15.0"
+version = "1.16.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
+checksum = "91622ff5e7162018101f2fea40d6ebf4a78bbe5a49736a2020649edf9693679e"
 
 [[package]]
 name = "equivalent"
@@ -400,9 +400,9 @@ checksum = "8f42a60cbdf9a97f5d2305f08a87dc4e09308d1276d28c869c684d7777685682"
 
 [[package]]
 name = "js-sys"
-version = "0.3.98"
+version = "0.3.99"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "67df7112613f8bfd9150013a0314e196f4800d3201ae742489d999db2f979f08"
+checksum = "142bc4740e452c1e57ade0cbc129f139c9093e354346f0872ef985f4f5cf5f11"
 dependencies = [
  "cfg-if",
  "futures-util",
@@ -436,15 +436,15 @@ checksum = "32a66949e030da00e8c7d4434b251670a91556f4144941d37452769c25d58a53"
 
 [[package]]
 name = "log"
-version = "0.4.29"
+version = "0.4.32"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5e5032e24019045c762d3c0f28f5b6b8bbf38563a65908389bf7978758920897"
+checksum = "953f07c43838f8e6f9758cab68bf5bed85465e7587ebe0b823f1bcd81978ad3a"
 
 [[package]]
 name = "memchr"
-version = "2.8.0"
+version = "2.8.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "f8ca58f447f06ed17d5fc4043ce1b10dd205e060fb3ce5b979b8ed8e59ff3f79"
+checksum = "6b947ae49db0d222b1dbc6b113ce7248a3fc3a6ca21b696717bfc000ba4484d8"
 
 [[package]]
 name = "num-traits"
@@ -560,9 +560,9 @@ dependencies = [
 
 [[package]]
 name = "qtty"
-version = "0.8.1"
+version = "0.8.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "61c018436452aeeb42be95367e44877b9237508908ee2d0f0ea6c13bcb969dea"
+checksum = "311988508597047d023edfb10cf4299b22a317be2b2c95745aaa4d489d9a9265"
 dependencies = [
  "qtty-core",
  "qtty-derive",
@@ -570,9 +570,9 @@ dependencies = [
 
 [[package]]
 name = "qtty-core"
-version = "0.8.1"
+version = "0.8.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1e4049b4883c30d10d9febc7a24088c4a5690e0f9fbdb63d91ab679352040744"
+checksum = "dc07e03ed9d81d4b4494b661635dc87f0a1a38cdce3e51136c39af6887716b02"
 dependencies = [
  "libm",
  "qtty-derive",
@@ -582,9 +582,9 @@ dependencies = [
 
 [[package]]
 name = "qtty-derive"
-version = "0.8.1"
+version = "0.8.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5c3ddfe91f467976a8bd7cab95750e3ffab4eb1b2875bb939fd9466e02790300"
+checksum = "7ee2574d42b017be334315c8cf4fdb241d728b093c5c7aba9bb1e7c5a3a708e8"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -783,9 +783,9 @@ dependencies = [
 
 [[package]]
 name = "serde_json"
-version = "1.0.149"
+version = "1.0.150"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "83fc039473c5595ace860d8c4fafa220ff474b3fc6bfdb4293327f1a37e94d86"
+checksum = "e8014e44b4736ed0538adeecded0fce2a272f22dc9578a7eb6b2d9993c74cfb9"
 dependencies = [
  "itoa",
  "memchr",
@@ -805,9 +805,9 @@ dependencies = [
 
 [[package]]
 name = "shlex"
-version = "1.3.0"
+version = "2.0.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64"
+checksum = "f8fadd59c855ef2080decdef8ff161eb6661b86933c9d82e5ba29dc602a55aba"
 
 [[package]]
 name = "slab"
@@ -920,9 +920,9 @@ dependencies = [
 
 [[package]]
 name = "typenum"
-version = "1.20.0"
+version = "1.20.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "40ce102ab67701b8526c123c1bab5cbe42d7040ccfd0f64af1a385808d2f43de"
+checksum = "b6f5e870be6c3b371b77fe0ee0bafb859fa4964b4404c27de1d380043c4dda20"
 
 [[package]]
 name = "unarray"
@@ -981,9 +981,9 @@ dependencies = [
 
 [[package]]
 name = "wasm-bindgen"
-version = "0.2.121"
+version = "0.2.122"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "49ace1d07c165b0864824eee619580c4689389afa9dc9ed3a4c75040d82e6790"
+checksum = "3ed04576f974d2b2fba0f38c51dbc5518011e38c36bf1143164be765528fd409"
 dependencies = [
  "cfg-if",
  "once_cell",
@@ -994,9 +994,9 @@ dependencies = [
 
 [[package]]
 name = "wasm-bindgen-macro"
-version = "0.2.121"
+version = "0.2.122"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "8e68e6f4afd367a562002c05637acb8578ff2dea1943df76afb9e83d177c8578"
+checksum = "916151b09da36bd82f6615cbf3a419e2f0ba23a03c6160e8e92eb6bd4aa1dec6"
 dependencies = [
  "quote",
  "wasm-bindgen-macro-support",
@@ -1004,9 +1004,9 @@ dependencies = [
 
 [[package]]
 name = "wasm-bindgen-macro-support"
-version = "0.2.121"
+version = "0.2.122"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d95a9ec35c64b2a7cb35d3fead40c4238d0940c86d107136999567a4703259f2"
+checksum = "299047362ccbfce148b67ab7e73349f77748e00c8296f9542adfad2ad82c5c5e"
 dependencies = [
  "bumpalo",
  "proc-macro2",
@@ -1017,9 +1017,9 @@ dependencies = [
 
 [[package]]
 name = "wasm-bindgen-shared"
-version = "0.2.121"
+version = "0.2.122"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c4e0100b01e9f0d03189a92b96772a1fb998639d981193d7dbab487302513441"
+checksum = "9a929b2c61f11ba3e9bc35b50c1f25cb38e0e892c0c231ae2b8cf78d5dad4437"
 dependencies = [
  "unicode-ident",
 ]
@@ -1060,9 +1060,9 @@ dependencies = [
 
 [[package]]
 name = "web-sys"
-version = "0.3.98"
+version = "0.3.99"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "4b572dff8bcf38bad0fa19729c89bb5748b2b9b1d8be70cf90df697e3a8f32aa"
+checksum = "6d621441cfc37b84979402712047321980c178f299193a3589d05b99e8763436"
 dependencies = [
  "js-sys",
  "wasm-bindgen",
@@ -1216,18 +1216,18 @@ dependencies = [
 
 [[package]]
 name = "zerocopy"
-version = "0.8.48"
+version = "0.8.50"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "eed437bf9d6692032087e337407a86f04cd8d6a16a37199ed57949d415bd68e9"
+checksum = "3b065d4f0e55f82fae73202e189638116a87c55ab6b8e6c2721e13dd9d854ad1"
 dependencies = [
  "zerocopy-derive",
 ]
 
 [[package]]
 name = "zerocopy-derive"
-version = "0.8.48"
+version = "0.8.50"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "70e3cd084b1788766f53af483dd21f93881ff30d7320490ec3ef7526d203bad4"
+checksum = "0b631b19d36a892ab55420c92dbc83ccd79274f25be714855d3074aa71cab639"
 dependencies = [
  "proc-macro2",
  "quote",

From 8773a20da3dc146f4c4e4340e4b92f42e42c97c0 Mon Sep 17 00:00:00 2001
From: VPRamon <vallespuigramon@gmail.com>
Date: Sun, 7 Jun 2026 14:59:29 +0200
Subject: [PATCH 5/5] test

---
 tests/roots.rs | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/tests/roots.rs b/tests/roots.rs
index d9f3346..fdbf0fc 100644
--- a/tests/roots.rs
+++ b/tests/roots.rs
@@ -205,8 +205,8 @@ fn fit_degree_one_produces_two_coefficients() {
 }
 
 #[test]
-fn chebyshev_t24_roots() {
-    let n = 24usize;
+fn chebyshev_t16_high_degree_roots() {
+    let n = 16usize;
     let mut coeffs = vec![0.0_f64; n + 1];
     coeffs[n] = 1.0;
     let p = ChebySeriesDyn::new(coeffs).unwrap();
@@ -216,6 +216,22 @@ fn chebyshev_t24_roots() {
     assert_residuals(&p, &roots, options.zero_tol);
 }
 
+#[test]
+fn high_degree_t24_may_miss_roots_but_residuals_hold() {
+    let n = 24usize;
+    let mut coeffs = vec![0.0_f64; n + 1];
+    coeffs[n] = 1.0;
+    let p = ChebySeriesDyn::new(coeffs).unwrap();
+    let options = opts();
+    let roots = p.roots_with(options);
+    assert!(
+        roots.len() >= n / 2,
+        "expected most T_24 roots, got {roots:?}"
+    );
+    assert!(roots.len() <= n);
+    assert_residuals(&p, &roots, options.zero_tol);
+}
+
 #[test]
 fn strict_zero_tol_can_exclude_approximate_roots() {
     let p = fit_dyn_from_fn(12, |x| (x - 0.3).powi(2)).unwrap();