fix conversion and root finding errors for chebyshev

Author: mileslucas

src/common.jl

@@ -32,7 +32,11 @@ julia> fromroots(r)
 Polynomial(6 - 5*x + x^2)
 ```
 """
-fromroots(P::Type{<:AbstractPolynomial}, r::AbstractVector; var::SymbolLike = :x)
+function fromroots(P::Type{<:AbstractPolynomial}, roots::AbstractVector; var::SymbolLike = :x)
+    x = variable(P, var)
+    p = [x - r for r in roots]
+    return truncate!(reduce(*, p))
+end
 fromroots(r::AbstractVector{<:Number}; var::SymbolLike = :x) = fromroots(Polynomial, r, var = var)
 fromroots(r; var::SymbolLike = :x) = fromroots(collect(r), var = var)
 
@@ -100,12 +104,12 @@ function roots(p::AbstractPolynomial{T}) where {T <: Number}
     if d < 1 return [] end
     d == 1 && return [-p[0] / p[1]]
 
-    chopped_trimmed = chop(truncate(p))
+    chopped_trimmed = truncate(p)
     n_trail = length(p) - length(chopped_trimmed)
     comp = companion(chopped_trimmed)
     L = eigvals(rot180(comp))
-    append!(L, zeros(n_trail))
-    return sort!(L, rev = true, by = norm)
+    append!(L, zeros(eltype(L), n_trail))
+    return sort!(L, rev = true, by = real)
 end
 
 """
@@ -321,8 +325,8 @@ Given values of x that are assumed to be unbounded (-∞, ∞), return values re
 julia> x = -10:10
 -10:10
 
-julia> mapdomain(ChebyshevT, x)
--1.0:0.1:1.0
+julia> extrema(mapdomain(ChebyshevT, x))
+(-1.0, 1.0)
 
 ```
 """
@@ -383,8 +387,20 @@ identity
 =#
 Base.copy(p::P) where {P <: AbstractPolynomial} = P(copy(p.coeffs), p.var)
 Base.hash(p::AbstractPolynomial, h::UInt) = hash(p.var, hash(p.coeffs, h))
+"""
+    zero(::Type{<:AbstractPolynomial})
+    zero(::AbstractPolynomial)
+
+Returns a representation of 0 as the given polynomial.
+"""
 Base.zero(::Type{P}) where {P <: AbstractPolynomial} = P(zeros(1))
 Base.zero(p::P) where {P <: AbstractPolynomial} = zero(P)
+"""
+    one(::Type{<:AbstractPolynomial})
+    one(::AbstractPolynomial)
+
+Returns a representation of 1 as the given polynomial.
+"""
 Base.one(::Type{P}) where {P <: AbstractPolynomial} = P(ones(1))
 Base.one(p::P) where {P <: AbstractPolynomial} = one(P)
 
@@ -430,7 +446,7 @@ function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: Abstrac
 end
 
 """
-    gcd(a::Polynomial, b::Polynomial)
+    gcd(a::AbstractPolynomial, b::AbstractPolynomial)
 
 Find the greatest common denominator of two polynomials recursively using
 [Euclid's algorithm](http://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid.27s_algorithm).

src/polynomials/ChebyshevT.jl

@@ -3,7 +3,26 @@ export ChebyshevT
 """
     ChebyshevT{<:Number}(coeffs::AbstractVector, var=:x)
 
-Chebyshev polynomial of the first kind
+Chebyshev polynomial of the first kind.
+
+Construct a polynomial from its coefficients `a`, lowest order first, optionally in
+terms of the given variable `x`. `x` can be a character, symbol, or string.
+
+# Examples
+
+```jldoctest
+julia> c = ChebyshevT([1, 0, 3, 4])
+ChebyshevT([1, 0, 2, 4])
+
+julia> ChebyshevT([1, 2, 3], :s)
+ChebyshevT([1, 2, 3])
+
+julia> one(ChebyshevT)
+ChebyshevT(1.0)
+
+julia> convert(Polynomial, c)
+Polynomial(-2 - 12*x + 6*x^2 + 16*x^3)
+```
 """
 struct ChebyshevT{T <: Number} <: AbstractPolynomial{T}
     coeffs::Vector{T}
@@ -59,21 +78,6 @@ function (ch::ChebyshevT{T})(x::S) where {T,S}
     return R(c0 + c1 * x)
 end
 
-function fromroots(P::Type{<:ChebyshevT}, roots::AbstractVector{T}; var::SymbolLike = :x) where {T <: Number}
-    p = [P([-r, 1]) for r in roots]
-    n = length(p)
-    while n > 1
-        m, r = divrem(n, 2)
-        tmp = [p[i] * p[i + m] for i in 1:m]
-        if r > 0
-            tmp[1] *= p[end]
-        end
-        p = tmp
-        n = m
-    end
-    return truncate!(p[1])
-end
-
 function vander(P::Type{<:ChebyshevT}, x::AbstractVector{T}, n::Integer) where {T <: Number}
     A = Matrix{T}(undef, length(x), n + 1)
     A[:, 1] .= one(T)
@@ -136,14 +140,14 @@ function companion(p::ChebyshevT{T}) where T
     d == 1 && return diagm(0 => [-p[0] / p[1]])
     R = eltype(one(T) / one(T))
 
-    scl = append!([1.0], √5 .* ones(d - 1))
+    scl = vcat(1.0, √0.5 .* ones(R, d - 1))
 
-    diag = append!([√5], 0.5 .* ones(d - 2))
-    comp = diagm(-1 => diag,
-                  1 => diag)
+    diag = vcat(√0.5, ones(R, d - 2) ./ 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)
+    comp[:, end] .-= monics[1:d] .* scl ./ scl[end] ./ 2
+    return R.(comp)
 end
 
 function Base.:+(p1::ChebyshevT, p2::ChebyshevT)

test/ChebyshevT.jl

@@ -53,8 +53,15 @@ end
 end
 
 @testset "Roots" begin
-    @test fromroots(ChebyshevT, [-1, 0, 1]) == ChebyshevT([0, -0.25, 0, 0.25])
-    @test fromroots(ChebyshevT, [-1im, 1im]) ≈ ChebyshevT([1.5 + 0im, 0 + 0im, 0.5 + 0im])
+    r = [-1, 0, 1]
+    c = fromroots(ChebyshevT, r)
+    @test c == ChebyshevT([0, -0.25, 0, 0.25])
+    @test roots(c) ≈ sort(r, rev = true)
+
+    r = [-1im, 1im]
+    c = fromroots(ChebyshevT, r)
+    @test c ≈ ChebyshevT([1.5 + 0im, 0 + 0im, 0.5 + 0im])
+    @test roots(c) ≈ sort(r, by = real, rev = true)
 end
 
 @testset "Values" begin
@@ -67,7 +74,9 @@ end
 @testset "Conversions" begin
     c1 = ChebyshevT([2.5, 2.0, 1.5])
     @test c1 ≈ Polynomial([1, 2, 3])
-    @test convert(Polynomial{Float64}, c1) == Polynomial{Float64}([1, 2, 3])
+    p = convert(Polynomial{Float64}, c1)
+    @test p == Polynomial{Float64}([1, 2, 3])
+    @test convert(ChebyshevT, p) == c1
 
 end
 
@@ -165,7 +174,7 @@ end
         @test derivative(cint) == cheb
     end
 end
-@info ""
+
 @testset "z-series" for i in 0:5
     # c to z
     input = vcat(2, ones(i))