improve ImmutablePolynomial

Author: jverzani

src/abstract.jl

@@ -13,6 +13,9 @@ An abstract container for various polynomials.
 """
 abstract type AbstractPolynomial{T} end
 
+# We want  ⟒(P{α…,T}) = P{α…}; this default
+# works for most cases
+⟒(P::Type{<:AbstractPolynomial}) = constructorof(P)
 
 """
     Polynomials.@register(name)
@@ -36,9 +39,9 @@ macro register(name)
     quote
         Base.convert(::Type{P}, p::P) where {P<:$poly} = p
         Base.convert(P::Type{<:$poly}, p::$poly{T}) where {T} = P(coeffs(p), p.var)
-        Base.promote_rule(::Type{$poly{T}}, ::Type{$poly{S}}) where {T,S} =
+        Base.promote_rule(::Type{<:$poly{T}}, ::Type{<:$poly{S}}) where {T,S} =
             $poly{promote_type(T, S)}
-        Base.promote_rule(::Type{$poly{T}}, ::Type{S}) where {T,S<:Number} =
+        Base.promote_rule(::Type{<:$poly{T}}, ::Type{S}) where {T,S<:Number} =
             $poly{promote_type(T, S)}
         $poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} =
             $poly{T}(coeffs, Symbol(var))

src/common.jl

@@ -172,7 +172,11 @@ In-place version of [`truncate`](@ref)
 """
 function truncate!(p::AbstractPolynomial{T};
     rtol::Real = Base.rtoldefault(real(T)),
-    atol::Real = 0,) where {T}
+                   atol::Real = 0,) where {T}
+
+    # 
+
+    
     max_coeff = maximum(abs, coeffs(p))
     thresh = max_coeff * rtol + atol
     map!(c->abs(c) <= thresh ? zero(T) : c, coeffs(p), coeffs(p))
@@ -392,7 +396,7 @@ Base.broadcastable(p::AbstractPolynomial) = Ref(p)
 function basis(::Type{P}, k::Int, _var::SymbolLike=:x; var=_var) where {P <: AbstractPolynomial}
     zs = zeros(Int, k+1)
     zs[end] = 1
-    P(zs, var)
+    ⟒(P){eltype(P)}(zs, var)
 end
 basis(p::P, k::Int, _var::SymbolLike=:x; var=_var) where {P<:AbstractPolynomial} = basis(P, k, var)
 
@@ -452,15 +456,15 @@ Base.hash(p::AbstractPolynomial, h::UInt) = hash(p.var, hash(coeffs(p), h))
 
 Returns a representation of 0 as the given polynomial.
 """
-Base.zero(::Type{P}, var=:x) where {P <: AbstractPolynomial} = P(zeros(eltype(P), 1), var)
+Base.zero(::Type{P}, var=:x) where {P <: AbstractPolynomial} = ⟒(P)(zeros(eltype(P), 1), var)
 Base.zero(p::P) where {P <: AbstractPolynomial} = zero(P, p.var)
 """
     one(::Type{<:AbstractPolynomial})
     one(::AbstractPolynomial)
 
 Returns a representation of 1 as the given polynomial.
 """
-Base.one(::Type{P}, var=:x) where {P <: AbstractPolynomial} = P(ones(eltype(P),1), var)
+Base.one(::Type{P}, var=:x) where {P <: AbstractPolynomial} = ⟒(P)(ones(eltype(P),1), var)
 Base.one(p::P) where {P <: AbstractPolynomial} = one(P, p.var)
 
 Base.oneunit(::Type{P}, args...) where {P <: AbstractPolynomial} = one(P, args...)
@@ -556,11 +560,11 @@ Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) =
 Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && p[0] == n
 Base.:(==)(n::Number, p::AbstractPolynomial) = p == n
 
+
 function Base.isapprox(p1::AbstractPolynomial{T},
     p2::AbstractPolynomial{S};
     rtol::Real = (Base.rtoldefault(T, S, 0)),
-                       atol::Real = 0,) where {T,S}
-    
+    atol::Real = 0,) where {T,S}
     p1, p2 = promote(p1, p2)
     check_same_variable(p1, p2)  || error("p1 and p2 must have same var")
     p1t = truncate(p1; rtol = rtol, atol = atol)
@@ -569,7 +573,6 @@ function Base.isapprox(p1::AbstractPolynomial{T},
         return false
     end
     isapprox(coeffs(p1t), coeffs(p2t), rtol = rtol, atol = atol)
-    
 end
 
 function Base.isapprox(p1::AbstractPolynomial{T},
@@ -583,6 +586,7 @@ function Base.isapprox(p1::AbstractPolynomial{T},
     isapprox(coeffs(p1t), [n], rtol = rtol, atol = atol)
 end
 
+
 Base.isapprox(n::S,
     p1::AbstractPolynomial{T};
     rtol::Real = (Base.rtoldefault(T, S, 0)),

src/contrib.jl

@@ -33,11 +33,14 @@ end
 ## Slight modification when `x` is a matrix
 ## Remove once dependencies for Julia 1.0.0 are dropped
 function evalpoly(x::S, p::Tuple) where {S}
+    p == () && return zero(S)
     if @generated
         N = length(p.parameters)
         ex = :(p[end]*_one(S))
-        for i in N-1:-1:1
-            ex = :(_muladd(x, $ex, p[$i]))
+        if N > 0
+            for i in N-1:-1:1
+                ex = :(_muladd(x, $ex, p[$i]))
+            end
         end
         ex
     else