From 080650daa8d3d1afba5f869eb66d45c9c2ce0b00 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 09:25:44 -0500
Subject: [PATCH 01/16] roots of chebyshev polynomials

---
 Project.toml          |  4 +++-
 src/FastChebInterp.jl |  1 +
 src/roots.jl          | 42 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 46 insertions(+), 1 deletion(-)
 create mode 100644 src/roots.jl

diff --git a/Project.toml b/Project.toml
index 8ad34aa..a999fad 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,14 +3,16 @@ uuid = "cf66c380-9a80-432c-aff8-4f9c79c0bdde"
 version = "1.2.0"
 
 [deps]
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 
 [compat]
 ChainRulesCore = "1"
 ChainRulesTestUtils = "1"
 FFTW = "1.0"
+LinearAlgebra = "<0.0.1, 1"
 StaticArrays = "0.12, 1.0"
 julia = "1.3"
 
diff --git a/src/FastChebInterp.jl b/src/FastChebInterp.jl
index 59535fa..453dc1e 100644
--- a/src/FastChebInterp.jl
+++ b/src/FastChebInterp.jl
@@ -58,6 +58,7 @@ Base.zero(c::ChebPoly{N,T,Td}) where {N,T,Td} = ChebPoly{N,T,Td}(zero(c.coefs),
 include("interp.jl")
 include("regression.jl")
 include("eval.jl")
+include("roots.jl")
 include("chainrules.jl")
 
 end # module
diff --git a/src/roots.jl b/src/roots.jl
new file mode 100644
index 0000000..7d95637
--- /dev/null
+++ b/src/roots.jl
@@ -0,0 +1,42 @@
+# Computing roots of Chebyshev real-valued polynomials in 1d via colleague matrices on
+# recursively subdivided intervals, as described in Trefethen, "Approximation
+# Theory and Approximation Practice", chapter 18.  https://epubs.siam.org/doi/10.1137/1.9781611975949
+#
+# Some additional implementation tricks were inspired by the chebfun source
+# code: https://github.com/chebfun/chebfun/blob/7574c77680d7e82b79626300bf255498271a72df/%40chebtech/roots.m
+
+using LinearAlgebra
+export colleague_matrix
+
+# colleague matrix for 1d array of coefficients, assuming
+# length > 2 and trailing zero coefficients have already been dropped.
+function _colleague_matrix(coefs::AbstractVector{<:Real})
+    n = length(coefs)
+    n > 2 || throw(ArgumentError("length(coefs) = $(length(coefs)) must be > 2"))
+    iszero(coefs[end]) && throw(ArgumentError("trailing coefficient must be nonzero"))
+
+    # colleage matrix, transposed to make it upper-Hessenberg (which doesn't change eigenvalues)
+    halves = fill(one(eltype(coefs))/2, n-2)
+    C = diagm(1 => halves, -1 => halves)
+    C[2,1] = 1
+    @views C[:,end] .-= coefs[1:end-1] ./ 2coefs[end]
+    return C
+end
+
+"""
+    colleague_matrix(c::ChebPoly{1,<:Real}; tol=5eps)
+
+Return the "colleague matrix" whose eigenvalues are the roots of the 1d real-valued
+Chebyshev polynomial `c`, dropping trailing coefficients whose relative contributions
+are `< tol` (defaulting to 5 × floating-point `eps`).   Throws an error if the resulting
+polynomial degree is less than 2.
+"""
+function colleague_matrix(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs))
+    abstol = sum(infnorm, c.coefs) * tol # absolute tolerance = L1 norm * tol
+    n = something(findlast(c -> infnorm(c) ≥ abstol, c.coefs), 1)
+    C = _colleague_matrix(@view c.coefs[1:n])
+    # scale and shift from [-1,1] to [lb,ub] via C * (ub-lb)/2 + (ub+lb)/2 * I
+    C .*= (c.ub[1] - c.lb[1])/2
+    diagview(C) .+= (c.ub[1] + c.lb[1])/2
+    return C
+end

From ac7729ea18bf6c3a15ce26ad40bc0442f8846124 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 09:52:39 -0500
Subject: [PATCH 02/16] handle degree < 2 cases

---
 src/roots.jl | 25 ++++++++++++++-----------
 1 file changed, 14 insertions(+), 11 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index 7d95637..a0d3db9 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -8,19 +8,23 @@
 using LinearAlgebra
 export colleague_matrix
 
-# colleague matrix for 1d array of coefficients, assuming
-# length > 2 and trailing zero coefficients have already been dropped.
+# colleague matrix for 1d array of Chebyshev coefficients, assuming
+# trailing zero coefficients have already been dropped.
 function _colleague_matrix(coefs::AbstractVector{<:Real})
     n = length(coefs)
-    n > 2 || throw(ArgumentError("length(coefs) = $(length(coefs)) must be > 2"))
+    n <= 1 && return float(eltype(coefs))[;;] # 0×0 case (no roots)
     iszero(coefs[end]) && throw(ArgumentError("trailing coefficient must be nonzero"))
 
-    # colleage matrix, transposed to make it upper-Hessenberg (which doesn't change eigenvalues)
-    halves = fill(one(eltype(coefs))/2, n-2)
-    C = diagm(1 => halves, -1 => halves)
-    C[2,1] = 1
-    @views C[:,end] .-= coefs[1:end-1] ./ 2coefs[end]
-    return C
+    if n == 2 # trivial 1×1 (degree 1) case
+        return [float(-coefs[1])/coefs[2];;]
+    else
+        # colleague matrix, transposed to make it upper-Hessenberg (which doesn't change eigenvalues)
+        halves = fill(one(float(eltype(coefs)))/2, n-2)
+        C = diagm(1 => halves, -1 => halves)
+        C[2,1] = 1
+        @views C[:,end] .-= coefs[1:end-1] ./ 2coefs[end]
+        return C
+    end
 end
 
 """
@@ -28,8 +32,7 @@ end
 
 Return the "colleague matrix" whose eigenvalues are the roots of the 1d real-valued
 Chebyshev polynomial `c`, dropping trailing coefficients whose relative contributions
-are `< tol` (defaulting to 5 × floating-point `eps`).   Throws an error if the resulting
-polynomial degree is less than 2.
+are `< tol` (defaulting to 5 × floating-point `eps`).
 """
 function colleague_matrix(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs))
     abstol = sum(infnorm, c.coefs) * tol # absolute tolerance = L1 norm * tol

From da0d99ef1fae07e43adffc946908387f0d4e9345 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 12:43:56 -0500
Subject: [PATCH 03/16] added roots

---
 src/roots.jl | 64 +++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 53 insertions(+), 11 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index a0d3db9..3f64f9f 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -6,11 +6,11 @@
 # code: https://github.com/chebfun/chebfun/blob/7574c77680d7e82b79626300bf255498271a72df/%40chebtech/roots.m
 
 using LinearAlgebra
-export colleague_matrix
+export colleague_matrix, roots
 
 # colleague matrix for 1d array of Chebyshev coefficients, assuming
 # trailing zero coefficients have already been dropped.
-function _colleague_matrix(coefs::AbstractVector{<:Real})
+function _colleague_matrix(coefs::AbstractVector{<:Number})
     n = length(coefs)
     n <= 1 && return float(eltype(coefs))[;;] # 0×0 case (no roots)
     iszero(coefs[end]) && throw(ArgumentError("trailing coefficient must be nonzero"))
@@ -27,19 +27,61 @@ function _colleague_matrix(coefs::AbstractVector{<:Real})
     end
 end
 
+function _colleague_matrix(coefs::AbstractVector{<:Number}, lb::Real, ub::Real)
+    C = _colleague_matrix(coefs)
+    # scale and shift from [-1,1] to [lb,ub] via C * (ub-lb)/2 + (ub+lb)/2 * I
+    C .*= (ub - lb)/2
+    diagview(C) .+= (ub + lb)/2
+    return C
+end
+
 """
-    colleague_matrix(c::ChebPoly{1,<:Real}; tol=5eps)
+    colleague_matrix(c::ChebPoly{1,<:Number}; tol=5eps)
 
 Return the "colleague matrix" whose eigenvalues are the roots of the 1d real-valued
 Chebyshev polynomial `c`, dropping trailing coefficients whose relative contributions
 are `< tol` (defaulting to 5 × floating-point `eps`).
 """
-function colleague_matrix(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs))
-    abstol = sum(infnorm, c.coefs) * tol # absolute tolerance = L1 norm * tol
-    n = something(findlast(c -> infnorm(c) ≥ abstol, c.coefs), 1)
-    C = _colleague_matrix(@view c.coefs[1:n])
-    # scale and shift from [-1,1] to [lb,ub] via C * (ub-lb)/2 + (ub+lb)/2 * I
-    C .*= (c.ub[1] - c.lb[1])/2
-    diagview(C) .+= (c.ub[1] + c.lb[1])/2
-    return C
+function colleague_matrix(c::ChebPoly{1,<:Number}; tol::Real=5*epsvals(c.coefs))
+    abstol = sum(abs, c.coefs) * tol # absolute tolerance = L1 norm * tol
+    n = something(findlast(c -> abs(c) ≥ abstol, c.coefs), 1)
+    return UpperHessenberg(_colleague_matrix(@view(c.coefs[1:n]), c.lb[1], c.ub[1]))
+end
+
+function filter_roots(c::ChebPoly{1,<:Real}, roots::AbstractVector{<:Number})
+    @inbounds lb, ub = c.lb[1], c.ub[1]
+    htol = eps(float(ub - lb)) * 100 # similar to chebfun
+    return [ clamp(real(r), lb, ub) for r in roots if abs(imag(r)) < htol && lb - htol <= real(r) <= ub + htol ]
+end
+
+if isdefined(LinearAlgebra.LAPACK, :hseqr!)
+    # see also LinearAlgebra.jl#1557 - optimized in-place eigenvalues for upper-Hessenberg colleague matrix
+    function hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T<:Union{LinearAlgebra.BlasReal, LinearAlgebra.BlasComplex}}
+        ilo, ihi, _ = LinearAlgebra.LAPACK.gebal!('S', triu!(C.data, -1))
+        return sort!(LinearAlgebra.LAPACK.hseqr!('E', 'N', 1, size(C,1), C.data, C.data)[3], by=reim)
+    end
+end
+hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T} = eigvals!(C.data)
+
+function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Integer=50)
+    tol > 0 || throw(ArgumentError("tolerance $tol for truncating coefficients must be > 0"))
+    abstol = sum(abs, c.coefs) * tol # absolute tolerance = L1 norm * tol
+    n = something(findlast(c -> abs(c) ≥ abstol, c.coefs), 1)
+    if n <= maxsize
+        λ = hesseneigvals!(UpperHessenberg(_colleague_matrix(@view c.coefs[1:n])))
+        λ .= (λ .+ 1) .* ((c.ub[1] - c.lb[1])/2) .+ c.lb[1] # scale and shift to [lb,ub]
+        return filter_roots(c, λ)
+    else
+        # roughly halve the domain, constructing new Chebyshev polynomials on each half, and
+        # call roots recursively.  as for chebfun, we split at an arbitrary point 0.004849834917525
+        # on [-1,1] rather than at 0 to avoid introducing additional spurious roots (since 0 is
+        # often a special point by symmetry).
+        split = oftype(tol, 1.004849834917525) * ((c.ub[1] - c.lb[1])/2) + c.lb[1]
+
+        # pick a fast order for the recursive DCT, should be highly composite, say 2^m
+        order = nextpow(2, length(c.coefs)-1)
+        c1 = chebinterp(c, order, c.lb[1], split; tol=tol)
+        c2 = chebinterp(c, order, split, c.ub[1]; tol=tol)
+        return vcat(roots(c1; tol=2tol, maxsize=maxsize), roots(c2; tol=2tol, maxsize=maxsize))
+    end
 end

From a4270193187b537d4e4bd244d6056972b90086bc Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 12:44:34 -0500
Subject: [PATCH 04/16] tweak

---
 src/roots.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index 3f64f9f..ca89ac6 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -73,7 +73,7 @@ function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Int
         return filter_roots(c, λ)
     else
         # roughly halve the domain, constructing new Chebyshev polynomials on each half, and
-        # call roots recursively.  as for chebfun, we split at an arbitrary point 0.004849834917525
+        # call roots recursively.  Following chebfun, we split at an arbitrary point 0.004849834917525
         # on [-1,1] rather than at 0 to avoid introducing additional spurious roots (since 0 is
         # often a special point by symmetry).
         split = oftype(tol, 1.004849834917525) * ((c.ub[1] - c.lb[1])/2) + c.lb[1]

From 2e85ed15adc4909bbe268290115b7ebac881aef3 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 12:53:01 -0500
Subject: [PATCH 05/16] docs for roots

---
 README.md    |  3 +++
 src/roots.jl | 12 ++++++++++++
 2 files changed, 15 insertions(+)

diff --git a/README.md b/README.md
index adf3952..da35f61 100644
--- a/README.md
+++ b/README.md
@@ -31,6 +31,9 @@ The FastChebInterp package also supports complex and vector-valued functions `f`
 the latter case, `c(y)` returns a vector of interpolants, and one can use `chebjacobian(c,y)`
 to compute the tuple `(c(y), J(y))` of the interpolant and its Jacobian matrix at `y`.
 
+There is also a function `roots(c)` to return an array of the real roots of `c` on the interval `[lb, ub]`,
+computed efficiently by a ["colleague"-matrix method (combined with recursive domain decomposition)](https://epubs.siam.org/doi/10.1137/1.9781611975949.ch18).
+
 ### Regression from arbitrary points
 
 We also have a function `chebregression(x, y, [lb, ub,], order)` that
diff --git a/src/roots.jl b/src/roots.jl
index ca89ac6..4b67853 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -63,6 +63,18 @@ if isdefined(LinearAlgebra.LAPACK, :hseqr!)
 end
 hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T} = eigvals!(C.data)
 
+"""
+    roots(c::ChebPoly{1,<:Real}; tol=5eps, maxsize::Integer=50)
+
+Returns an array of the real roots of `c` on the interval `[lb,ub]` (the lower and
+upper bounds of the interpolation domain).
+
+Uses a colleague-matrix method combined with recursive subdivision of the domain
+to keep the maximum matrix size `≤ maxsize`, following an algorithm described by
+Trefethen (*Approximation Theory and Approximation Practice*, ch. 18).  The `tol`
+argument is a relative tolerance for dropping small polynomial coefficients, defaulting
+to `5eps` where `eps` is the precision of `c`.
+"""
 function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Integer=50)
     tol > 0 || throw(ArgumentError("tolerance $tol for truncating coefficients must be > 0"))
     abstol = sum(abs, c.coefs) * tol # absolute tolerance = L1 norm * tol

From 8f83390c504236aae5f07d2673357e75b20b9f85 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 12:59:27 -0500
Subject: [PATCH 06/16] added tests for roots

---
 src/roots.jl     | 2 +-
 test/runtests.jl | 5 +++++
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index 4b67853..a59b07e 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -66,7 +66,7 @@ hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T} = eigvals!(C.data)
 """
     roots(c::ChebPoly{1,<:Real}; tol=5eps, maxsize::Integer=50)
 
-Returns an array of the real roots of `c` on the interval `[lb,ub]` (the lower and
+Returns a sorted array of the real roots of `c` on the interval `[lb,ub]` (the lower and
 upper bounds of the interpolation domain).
 
 Uses a colleague-matrix method combined with recursive subdivision of the domain
diff --git a/test/runtests.jl b/test/runtests.jl
index 6f771d2..5ba8915 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,6 +24,11 @@ Random.seed!(314159) # make chainrules tests deterministic
         @test chebgradient(interp, x1) ≈′ (f(x1), f′(x1))
         test_frule(interp, x1, rtol=sqrt(eps(T)), atol=sqrt(eps(T)))
         test_rrule(interp, x1, rtol=sqrt(eps(T)), atol=sqrt(eps(T)))
+
+        @test roots(chebinterp(x -> 1.1, 10, 0,1)) ≈ Float64[]
+        @test roots(chebinterp(x -> x - 0.1, 10, 0,1)) ≈ [0.1]
+        @test roots(chebinterp(x -> (x - 0.1)*(x-0.2), 10, 0,1)) ≈ [0.1, 0.2]
+        @test roots(chebinterp(cos, 10000, 0, 1000))/pi ≈ (0:317) .+ T(0.5)
     end
 end
 

From 5308cb1674fd40064b3559b8ce6976c02bf28605 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:04:23 -0500
Subject: [PATCH 07/16] document and test colleague_matrix

---
 README.md        | 3 ++-
 test/runtests.jl | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/README.md b/README.md
index da35f61..f59fefb 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,8 @@ the latter case, `c(y)` returns a vector of interpolants, and one can use `chebj
 to compute the tuple `(c(y), J(y))` of the interpolant and its Jacobian matrix at `y`.
 
 There is also a function `roots(c)` to return an array of the real roots of `c` on the interval `[lb, ub]`,
-computed efficiently by a ["colleague"-matrix method (combined with recursive domain decomposition)](https://epubs.siam.org/doi/10.1137/1.9781611975949.ch18).
+computed efficiently by a ["colleague"-matrix method (combined with recursive domain decomposition)](https://epubs.siam.org/doi/10.1137/1.9781611975949.ch18).  One can also obtain the colleague matrix (whose eigenvalues are the roots)
+directly via `colleague_matrix(c)` (which can be useful to examine complex roots and roots slightly outside `[lb, ub]`).
 
 ### Regression from arbitrary points
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 5ba8915..e85e531 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -29,6 +29,7 @@ Random.seed!(314159) # make chainrules tests deterministic
         @test roots(chebinterp(x -> x - 0.1, 10, 0,1)) ≈ [0.1]
         @test roots(chebinterp(x -> (x - 0.1)*(x-0.2), 10, 0,1)) ≈ [0.1, 0.2]
         @test roots(chebinterp(cos, 10000, 0, 1000))/pi ≈ (0:317) .+ T(0.5)
+        @test colleague_matrix(chebinterp(x -> (x - 0.1)*(x-0.2), 10, 0,1)) ≈ T[0.5 -0.24; 0.5 -0.2]
     end
 end
 

From 407825e6590e673094b9c418fb27cab9584fe686 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:06:34 -0500
Subject: [PATCH 08/16] fix empty-matrix syntax for Julia 1.3

---
 src/roots.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index a59b07e..57a4150 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -12,7 +12,7 @@ export colleague_matrix, roots
 # trailing zero coefficients have already been dropped.
 function _colleague_matrix(coefs::AbstractVector{<:Number})
     n = length(coefs)
-    n <= 1 && return float(eltype(coefs))[;;] # 0×0 case (no roots)
+    n <= 1 && return zeros(float(eltype(coefs)),0,0) # 0×0 case (no roots)
     iszero(coefs[end]) && throw(ArgumentError("trailing coefficient must be nonzero"))
 
     if n == 2 # trivial 1×1 (degree 1) case

From 051bd19ba39a9e2d5dc5af91c2dcfcf6aab63f6a Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:09:38 -0500
Subject: [PATCH 09/16] fix matrix-literal syntax for Julia 1.3

---
 src/roots.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index 57a4150..4999d3e 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -12,11 +12,13 @@ export colleague_matrix, roots
 # trailing zero coefficients have already been dropped.
 function _colleague_matrix(coefs::AbstractVector{<:Number})
     n = length(coefs)
-    n <= 1 && return zeros(float(eltype(coefs)),0,0) # 0×0 case (no roots)
+    n <= 1 && return Matrix{float(eltype(coefs))}(undef,0,0) # 0×0 case (no roots)
     iszero(coefs[end]) && throw(ArgumentError("trailing coefficient must be nonzero"))
 
     if n == 2 # trivial 1×1 (degree 1) case
-        return [float(-coefs[1])/coefs[2];;]
+        C = Matrix{float(eltype(coefs))}(undef,1,1)
+        C[1,1] = float(-coefs[1])/coefs[2]
+        return C
     else
         # colleague matrix, transposed to make it upper-Hessenberg (which doesn't change eigenvalues)
         halves = fill(one(float(eltype(coefs)))/2, n-2)

From 0bfd17f0a4d078ef2df1c7331351f3a8ed788c2e Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:15:16 -0500
Subject: [PATCH 10/16] compat for diagview

---
 src/roots.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index 4999d3e..deda48d 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -29,6 +29,10 @@ function _colleague_matrix(coefs::AbstractVector{<:Number})
     end
 end
 
+if !isdefined(Base, :diagview) # added in Julia 1.12
+    diagview(A::AbstractMatrix, k::Integer=0) = @view A[diagind(A, k, IndexStyle(A))]
+end
+
 function _colleague_matrix(coefs::AbstractVector{<:Number}, lb::Real, ub::Real)
     C = _colleague_matrix(coefs)
     # scale and shift from [-1,1] to [lb,ub] via C * (ub-lb)/2 + (ub+lb)/2 * I
@@ -56,7 +60,7 @@ function filter_roots(c::ChebPoly{1,<:Real}, roots::AbstractVector{<:Number})
     return [ clamp(real(r), lb, ub) for r in roots if abs(imag(r)) < htol && lb - htol <= real(r) <= ub + htol ]
 end
 
-if isdefined(LinearAlgebra.LAPACK, :hseqr!)
+if isdefined(LinearAlgebra.LAPACK, :hseqr!) # added in Julia 1.10
     # see also LinearAlgebra.jl#1557 - optimized in-place eigenvalues for upper-Hessenberg colleague matrix
     function hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T<:Union{LinearAlgebra.BlasReal, LinearAlgebra.BlasComplex}}
         ilo, ihi, _ = LinearAlgebra.LAPACK.gebal!('S', triu!(C.data, -1))

From 3e6277a65cf7c17021d6bd7d173898a8dabb87dd Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:23:01 -0500
Subject: [PATCH 11/16] grr, diagind also not defined in Julia 1.3

---
 src/roots.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index deda48d..83527b4 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -29,8 +29,12 @@ function _colleague_matrix(coefs::AbstractVector{<:Number})
     end
 end
 
-if !isdefined(Base, :diagview) # added in Julia 1.12
-    diagview(A::AbstractMatrix, k::Integer=0) = @view A[diagind(A, k, IndexStyle(A))]
+
+if !isdefined(LinearAlgebra, :diagview) # added in Julia 1.12
+    if !isdefined(LinearAlgebra, :diagind)
+        diagind(A::Matrix) = range(1, step=size(A,1)+1, length=min(size(A)...))
+    end
+    diagview(A::AbstracMatrixMatrix) = @view A[diagind(A)]
 end
 
 function _colleague_matrix(coefs::AbstractVector{<:Number}, lb::Real, ub::Real)

From 77492539dd03bd93450cbb64fb3a1d2b4f75fd29 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:23:19 -0500
Subject: [PATCH 12/16] whoops

---
 src/roots.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index 83527b4..2682bfe 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -34,7 +34,7 @@ if !isdefined(LinearAlgebra, :diagview) # added in Julia 1.12
     if !isdefined(LinearAlgebra, :diagind)
         diagind(A::Matrix) = range(1, step=size(A,1)+1, length=min(size(A)...))
     end
-    diagview(A::AbstracMatrixMatrix) = @view A[diagind(A)]
+    diagview(A::Matrix) = @view A[diagind(A)]
 end
 
 function _colleague_matrix(coefs::AbstractVector{<:Number}, lb::Real, ub::Real)

From 06def561648b7eac6afbb0d054cae776d8634451 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:26:05 -0500
Subject: [PATCH 13/16] paranoia in fallback hesseneigvals

---
 src/roots.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index 2682bfe..5ec0f0c 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -71,7 +71,7 @@ if isdefined(LinearAlgebra.LAPACK, :hseqr!) # added in Julia 1.10
         return sort!(LinearAlgebra.LAPACK.hseqr!('E', 'N', 1, size(C,1), C.data, C.data)[3], by=reim)
     end
 end
-hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T} = eigvals!(C.data)
+hesseneigvals!(C::UpperHessenberg{T,Matrix{T}}) where {T} = eigvals!(triu!(C.data, -1))
 
 """
     roots(c::ChebPoly{1,<:Real}; tol=5eps, maxsize::Integer=50)

From 1b32234cd3da959faa5d765e739b0741ad7f8251 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:28:27 -0500
Subject: [PATCH 14/16] tweak

---
 src/roots.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index 5ec0f0c..5b639a2 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -91,14 +91,14 @@ function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Int
     n = something(findlast(c -> abs(c) ≥ abstol, c.coefs), 1)
     if n <= maxsize
         λ = hesseneigvals!(UpperHessenberg(_colleague_matrix(@view c.coefs[1:n])))
-        λ .= (λ .+ 1) .* ((c.ub[1] - c.lb[1])/2) .+ c.lb[1] # scale and shift to [lb,ub]
+        @inbounds λ .= (λ .+ 1) .* ((c.ub[1] - c.lb[1])/2) .+ c.lb[1] # scale and shift to [lb,ub]
         return filter_roots(c, λ)
     else
         # roughly halve the domain, constructing new Chebyshev polynomials on each half, and
         # call roots recursively.  Following chebfun, we split at an arbitrary point 0.004849834917525
         # on [-1,1] rather than at 0 to avoid introducing additional spurious roots (since 0 is
         # often a special point by symmetry).
-        split = oftype(tol, 1.004849834917525) * ((c.ub[1] - c.lb[1])/2) + c.lb[1]
+        @inbounds split = oftype(float(c.lb[1]), 1.004849834917525) * ((c.ub[1] - c.lb[1])/2) + c.lb[1]
 
         # pick a fast order for the recursive DCT, should be highly composite, say 2^m
         order = nextpow(2, length(c.coefs)-1)

From 5bf19d61e70390696fb145339f145cd55f82ff3b Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:30:12 -0500
Subject: [PATCH 15/16] tweak

---
 src/roots.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/roots.jl b/src/roots.jl
index 5b639a2..04a6a41 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -101,7 +101,7 @@ function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Int
         @inbounds split = oftype(float(c.lb[1]), 1.004849834917525) * ((c.ub[1] - c.lb[1])/2) + c.lb[1]
 
         # pick a fast order for the recursive DCT, should be highly composite, say 2^m
-        order = nextpow(2, length(c.coefs)-1)
+        order = nextpow(2, n-1)
         c1 = chebinterp(c, order, c.lb[1], split; tol=tol)
         c2 = chebinterp(c, order, split, c.ub[1]; tol=tol)
         return vcat(roots(c1; tol=2tol, maxsize=maxsize), roots(c2; tol=2tol, maxsize=maxsize))

From c6cc9552566903f70bcd918b1030e4c87331a674 Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@alum.mit.edu>
Date: Sat, 28 Feb 2026 13:32:06 -0500
Subject: [PATCH 16/16] consistent truncation

---
 src/roots.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index 04a6a41..7bfbc8b 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -102,8 +102,8 @@ function roots(c::ChebPoly{1,<:Real}; tol::Real=5*epsvals(c.coefs), maxsize::Int
 
         # pick a fast order for the recursive DCT, should be highly composite, say 2^m
         order = nextpow(2, n-1)
-        c1 = chebinterp(c, order, c.lb[1], split; tol=tol)
-        c2 = chebinterp(c, order, split, c.ub[1]; tol=tol)
+        c1 = chebinterp(c, order, c.lb[1], split; tol=2tol)
+        c2 = chebinterp(c, order, split, c.ub[1]; tol=2tol)
         return vcat(roots(c1; tol=2tol, maxsize=maxsize), roots(c2; tol=2tol, maxsize=maxsize))
     end
 end