companion working for ChebyshevT

Author: mileslucas

src/polynomials/ChebyshevT.jl

@@ -116,19 +116,21 @@ end
 # end
 
 ##
-# function companion(p::ChebyshevT{T}) where T
-#     d = length(p) - 1
-#     d < 1 && error("Series must have degree greater than 1")
-#     d == 1 && return diagm(0 => [-p[0] / p[1]])
-#     R = eltype(one(T) / p.coeffs[end])
-
-#     comp = diagm(-1 => ones(d - 1),
-#                   1 => ones(d - 1))
-#     comp[2, 1] = 2
-#     monics = p.coeffs ./ p.coeffs[end]
-#     comp[:, end] .= -monics[1:d]
-#     return R.(comp ./ 2)
-# end
+function companion(p::ChebyshevT{T}) where T
+    d = length(p) - 1
+    d < 1 && error("Series must have degree greater than 1")
+    d == 1 && return diagm(0 => [-p[0] / p[1]])
+    R = eltype(one(T) / one(T))
+
+    scl = append!([1.0], √5 .* ones(d - 1))
+
+    diag = append!([√5], 0.5 .* ones(d - 2))
+    comp = diagm(-1 => diag,
+                  1 => diag)
+    monics = p.coeffs ./ p.coeffs[end]
+    comp[:, end] .-= monics[1:d] .* scl ./ scl[end]
+    return R.(comp ./ 2)
+end
 
 function Base.:+(p1::ChebyshevT, p2::ChebyshevT)
     p1.var != p2.var && error("Polynomials must have same variable")

test/ChebyshevT.jl

@@ -71,6 +71,20 @@ end
 
 end
 
+@testset "Companion" begin
+    c_null = ChebyshevT(Int[])
+    c_1 = ChebyshevT([1])
+    @test_throws ErrorException companion(c_null)
+    @test_throws ErrorException companion(c_1)
+    for i in 1:5
+        coef = push!(zeros(i), 1)
+        c = ChebyshevT(coef)
+        @test size(companion(c)) == (i, i)
+    end
+    c = ChebyshevT([1, 2])
+    @test companion(c)[1, 1] == -0.5
+end
+
 @testset "Arithmetic $i, $j" for i in 1:5, j in 1:5
     # multiplication
     target = zeros(i + j + 1)