Here is a quick implementation based on the paper by Amiraslami, Corless and Lancaster

```
using LinearAlgebra, SpecialPolynomials
function generalized_comrade(coeffs::Vector{T},OP) where T <: Number
p = length(coeffs)
an_p = SpecialPolynomials.an.(OP,collect(0:p-2))
bn_p = SpecialPolynomials.bn.(OP,collect(0:p-1))
cn_p = SpecialPolynomials.cn.(OP,collect(1:p-1))
kn_p = SpecialPolynomials.leading_term.(OP,[p-1,p])
B = Diagonal([ones(T,p-2);kn_p[2]*coeffs[end]])
A = zeros(T,p-1,p-1)
A[1,1] = bn_p[1]
A[2,1] = an_p[1]
@inbounds for r in 2:p-2
A[r-1,r] = cn_p[r-1]
A[r ,r] = bn_p[r]
A[r+1,r] = an_p[r]
end
@views A[:,end] .= .-kn_p[1] * coeffs[1:end-1]
A[end-1,end] += kn_p[2] * cn_p[end] * coeffs[end]
A[end ,end] += kn_p[2] * bn_p[end] * coeffs[end]
return UpperHessenberg(A),B
end
```

There is however a non-negligible difference between this implementation and what is found by back substitution and companion matrix computation (which is what is currently done by SpecialPolynomials AFAICT).

I canâ€™t figure if this comes from a bug on my part (I checked again,but I may have missed something) or from the back substitution.

See :

```
using Polynomials
coeffs = rand(6)
L = Legendre(coeffs)
P = convert(Polynomial,L)
ref_roots = roots(P)
A,B = generalized_comrade(coeffs,Legendre)
comp_roots = eigvals(Matrix(A),Matrix(B)) # Don't know if diagonal and hessenberg structure can be leveraged there)
```

EDIT, I checked for Chebyshev polynomials, and the results check out w.r.t the analystical expression of the Chebyshev nodes

```
using Plots
n = 30
Cnodes(n) = cos.((2 .* (1:n) .- 1)./(2 .* n) .* pi)
C = basis(Chebyshev,n)
r = roots(C)
A,B = generalized_comrade(C.coeffs,Chebyshev)
rr = eigvals(Matrix(A),Matrix(B))
rel(a,b) = abs(a-b)/abs(b)
plot(rel.(-r,Cnodes(n)))
plot!(rel.(-rr,Cnodes(n)))
```