diff --git a/Project.toml b/Project.toml
index bb65eb73..7af3aee1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "2.0.12"
+version = "2.0.13"
 
 [deps]
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
diff --git a/docs/src/index.md b/docs/src/index.md
index 835e7c01..47616800 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -161,7 +161,7 @@ By design, this is not type-stable; the returned roots may be real or complex.
 
 The default `roots` function uses the eigenvalues of the
 [companion](https://en.wikipedia.org/wiki/Companion_matrix) matrix for
-a polynomial. This is an `𝑶(n^3)` operation. 
+a polynomial. This is an `𝑶(n^3)` operation.
 
 For polynomials with `BigFloat` coefficients, the
 `GenericLinearAlgebra` package can be seamlessly used:
@@ -204,7 +204,7 @@ julia> PolynomialRoots.roots(coeffs(p))
  1.0000000000000002 + 0.0im
 ```
 
-The roots are always returned as complex numbers. 
+The roots are always returned as complex numbers.
 
 
 * The
@@ -226,7 +226,7 @@ julia> AMRVW.roots(float.(coeffs(p)))
  2.9999999999999964 + 0.0im
 ```
 
-The roots are returned as complex numbers. 
+The roots are returned as complex numbers.
 
 Both `PolynomialRoots` and `AMRVW` are generic and work with
 `BigFloat` coefficients, for example.
@@ -429,7 +429,7 @@ Fit a polynomial (of degree `deg`) to `x` and `y` using polynomial interpolation
 using Plots, Polynomials
 xs = range(0, 10, length=10)
 ys = @. exp(-xs)
-f = fit(xs, ys) # degree = length(xs) - 1 
+f = fit(xs, ys) # degree = length(xs) - 1
 f2 = fit(xs, ys, 2) # degree = 2
 
 scatter(xs, ys, markerstrokewidth=0, label="Data")
@@ -583,7 +583,7 @@ Polynomial(4 + 2*x + 3*x^2)
 
 If `q` is non-constant, such as `variable(Polynomial, :y)`, then there would be an error due to the mismatched symbols. (The mathematical result would need a multivariate polynomial, not a univariate polynomial, as this package provides.)
 
-The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match. 
+The same conversion is done for polynomial multiplication: constant polynomials are treated as numbers; non-constant polynomials must have their symbols match.
 
 There is an oddity -- though the following two computations look the same, they are technically different:
 
@@ -595,7 +595,7 @@ julia> one(Polynomial, :y) + one(Polynomial, :x)
 Polynomial(2.0)
 ```
 
-Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`. 
+Both are constant polynomials over `Int`, but the first has the indeterminate `:y`, the second `:x`.
 
 This technical difference causes no issues with polynomial addition or multiplication, as there constant polynomials are treated as numbers, but can be an issue when constant polynomials are used as array elements.
 
@@ -669,6 +669,91 @@ julia> [1 p; p 1] + [1 2one(q); 3 4] # array{P{T,:x}} + array{P{T,:y}}
 Though were a non-constant polynomial with indeterminate `y` replacing
 `2one(q)` above, that addition would throw an error.
 
+
+## Non-number types for `T`
+
+The coefficients of the polynomial may be non-number types, such as matrices or other polynomials, albeit not every operation is fully supported.
+
+For example, a polynomial with matrix coefficients, might be constructed with:
+
+```jldoctest non_number
+julia> using Polynomials
+
+julia> a,b,c = [1 0;2 1], [1 0; 3 1], [1 0; 4 1]
+([1 0; 2 1], [1 0; 3 1], [1 0; 4 1])
+
+julia> p = Polynomial([a,b,c])
+Polynomial([1 0; 2 1] + [1 0; 3 1]*x + [1 0; 4 1]*x^2)
+
+julia> q = derivative(p)
+Polynomial([1 0; 3 1] + [2 0; 8 2]*x)
+```
+
+Various operations are available, `derivative` was shown above, here are the vector-space operations:
+
+```jldoctest non_number
+julia> 2p
+Polynomial([2 0; 4 2] + [2 0; 6 2]*x + [2 0; 8 2]*x^2)
+
+julia> p + q
+Polynomial([2 0; 5 2] + [3 0; 11 3]*x + [1 0; 4 1]*x^2)
+```
+
+polynomial multiplication:
+
+```jldoctest non_number
+julia> p * q
+Polynomial([1 0; 5 1] + [3 0; 18 3]*x + [3 0; 21 3]*x^2 + [2 0; 16 2]*x^3)
+```
+
+polynomial evaluation, here either with a scalar or a matrix:
+
+```jldoctest non_number
+julia> p(2)
+2×2 Matrix{Int64}:
+  7  0
+ 24  7
+
+julia> p(b)
+2×2 Matrix{Int64}:
+  3  0
+ 18  3
+```
+
+But if the type `T` lacks support of some generic functions, such as `zero(T)` and `one(T)`, then there may be issues. For example,  when `T <: AbstractMatrix` the output of `p-p` is an error, as the implementation assumes `zero(T)` is defined. For static arrays, this isn't an issue, as there is support for `zero(T)`. Other polynomial types, such as `SparsePolynomial` have less support, as some specialized methods assume more of the generic interface be implemented.
+
+Similarly, using polynomials for `T` is a possibility:
+
+```jldoctest non_number
+julia> a,b,c = Polynomial([1],:y), Polynomial([0,1],:y), Polynomial([0,0,1],:y)
+(Polynomial(1), Polynomial(y), Polynomial(y^2))
+
+julia> p = Polynomial([a,b,c], :x)
+Polynomial(Polynomial(1) + Polynomial(y)*x + Polynomial(y^2)*x^2)
+
+julia> q = derivative(p)
+Polynomial(Polynomial(y) + Polynomial(2*y^2)*x)
+```
+
+Again, much works:
+
+```jldoctest non_number
+julia> 2p
+Polynomial(Polynomial(2) + Polynomial(2*y)*x + Polynomial(2*y^2)*x^2)
+
+julia> p + q
+Polynomial(Polynomial(1 + y) + Polynomial(y + 2*y^2)*x + Polynomial(y^2)*x^2)
+
+julia> p(2)
+Polynomial(1 + 2*y + 4*y^2)
+
+julia> p(b)
+Polynomial(1 + y^2 + y^4)
+```
+
+But much doesn't. For example, implicit promotion can fail. For example, the scalar multiplication `p * b` will fail,  as the methods assume this is the fallback polynomial multiplication and not the intended scalar multiplication.
+
+
 ## Rational functions
 
 The package provides support for rational functions -- fractions of polynomials (for most types). The construction of the basic type mirrors the construction of rational numbers.
diff --git a/src/abstract.jl b/src/abstract.jl
index 1cd2aec2..a32c8ce4 100644
--- a/src/abstract.jl
+++ b/src/abstract.jl
@@ -1,17 +1,50 @@
 export AbstractPolynomial
 
+
+
 const SymbolLike = Union{AbstractString,Char,Symbol}
 
 """
     AbstractPolynomial{T,X}
 
-An abstract container for various polynomials. 
+An abstract type for various polynomials.
+
+A polynomial type holds an indeterminate `X`; coefficients of type `T`, stored in some container type; and an implicit basis, such as the standard basis.
 
 # Properties
 - `coeffs` - The coefficients of the polynomial
+
+# The type `T`
+
+`T` need not be `T <: Number`, at the moment it is not restricted
+
+Some `T`s will not be successful
+
+* scalar mult: `c::Number * p::Polynomial` should be defined
+* scalar mult: `c::T * p:Polynomial{T}` An  ambiguity when `T <: AbstractPolynomial`
+* scalar mult: `p:Polynomial{T} * c::T` need not commute
+
+* scalar add/sub: `p::Polynomial{T} + q::Polynomial{T}` should be defined
+* scalar sub: `p::Polynomial{T} - p::Polynomial{T}`  generally needs `zeros(T,1)` defined for `zero(Polynomial{T})`
+
+* poly mult: `p::Polynomial{T} * q::Polynomial{T}` Needs "`T * T`" defined (e.g. `Base.promote_op(*, Vector{Int}, Vector{Int}))` needs to be something.)
+* poly powers: `p::Polynomial{T}^2` needs "`T^2`" defined
+
+* implicit promotion: `p::Polynomial{T} + c::Number`  needs `convert(T, c)` defined
+* implicit promotion: `p::Polynomial{T} + c::T`  ambiguity if `T <: AbstractPolynomial`
+
+* evaluation: `p::Polynomial{T}(s::Number)`
+* evaluation `p::Polynomial{T}(c::T)`   needs `T*T` defined
+
+* derivatives: `derivative(p::Polynomial{T})`
+* integrals: `integrate(p::Polynomial{T})`
+
+
 """
 abstract type AbstractPolynomial{T,X} end
 
+
+## -----
 # We want  ⟒(P{α…,T}) = P{α…}; this default
 # works for most cases
 ⟒(P::Type{<:AbstractPolynomial}) = constructorof(P)
@@ -34,7 +67,7 @@ Polynomials.@register MyPolynomial
 ```
 
 # Implementations
-This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`). 
+This will implement simple self-conversions like `convert(::Type{MyPoly}, p::MyPoly) = p` and creates two promote rules. The first allows promotion between two types (e.g. `promote(Polynomial, ChebyshevT)`) and the second allows promotion between parametrized types (e.g. `promote(Polynomial{T}, Polynomial{S})`).
 
 For constructors, it implements the shortcut for `MyPoly(...) = MyPoly{T}(...)`, singleton constructor `MyPoly(x::Number, ...)`,  conversion constructor `MyPoly{T}(n::S, ...)`, and `variable` alternative  `MyPoly(var=:x)`.
 """
@@ -44,7 +77,7 @@ macro register(name)
         Base.convert(::Type{P}, p::P) where {P<:$poly} = p
         function Base.convert(P::Type{<:$poly}, p::$poly{T,X}) where {T,X}
             isconstant(p) && return constructorof(P){eltype(P),indeterminate(P)}(constantterm(p))
-            constructorof(P){eltype(P), indeterminate(P,p)}(coeffs(p))
+            constructorof(P){eltype(P), indeterminate(P,p)}(_coeffs(p))
         end
         Base.promote(p::P, q::Q) where {X, T, P <:$poly{T,X}, Q <: $poly{T,X}} = p,q
         Base.promote_rule(::Type{<:$poly{T,X}}, ::Type{<:$poly{S,X}}) where {T,S,X} =  $poly{promote_type(T, S),X}
@@ -52,7 +85,7 @@ macro register(name)
             $poly{promote_type(T, S),X}
         $poly(coeffs::AbstractVector{T}, var::SymbolLike = :x) where {T} =
             $poly{T, Symbol(var)}(coeffs)
-        $poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S<:Number} =
+        $poly{T}(x::AbstractVector{S}, var::SymbolLike = :x) where {T,S} =
             $poly{T,Symbol(var)}(T.(x))
         function $poly(coeffs::G, var::SymbolLike=:x) where {G}
             !Base.isiterable(G) && throw(ArgumentError("coeffs is not iterable"))
@@ -76,11 +109,11 @@ macro registerN(name, params...)
     αs = tuple(esc.(params)...)
     quote
         Base.convert(::Type{P}, q::Q) where {$(αs...),T, P<:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = q
-        Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q        
+        Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q
         Base.promote(p::P, q::Q) where {$(αs...),T, X, P <:$poly{$(αs...),T,X}, Q <: $poly{$(αs...),T,X}} = p,q
         Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{<:$poly{$(αs...),S,X}}) where {$(αs...),T,S,X} =
             $poly{$(αs...),promote_type(T, S),X}
-        Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{S}) where {$(αs...),T,X,S<:Number} = 
+        Base.promote_rule(::Type{<:$poly{$(αs...),T,X}}, ::Type{S}) where {$(αs...),T,X,S<:Number} =
             $poly{$(αs...),promote_type(T,S),X}
 
         function $poly{$(αs...),T}(x::AbstractVector{S}, var::SymbolLike = :x) where {$(αs...),T,S}
@@ -98,4 +131,3 @@ macro registerN(name, params...)
         (p::$poly)(x) = evalpoly(x, p)
     end
 end
-
diff --git a/src/common.jl b/src/common.jl
index db9d23bd..c473a88a 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -81,7 +81,7 @@ is [`Polynomial`](@ref).
 
 ## Weights
 
-Weights may be assigned to the points by specifying a vector or matrix of weights. 
+Weights may be assigned to the points by specifying a vector or matrix of weights.
 
 When specified as a vector, `[w₁,…,wₙ]`, the weights should be
 non-negative as the minimization problem is `argmin_β Σᵢ wᵢ |yᵢ - Σⱼ
@@ -89,9 +89,9 @@ Vᵢⱼ βⱼ|² = argmin_β || √(W)⋅(y - V(x)β)||²`, where, `W` the digon
 matrix formed from `[w₁,…,wₙ]`, is used for the solution, `V` being
 the Vandermonde matrix of `x` corresponding to the specified
 degree. This parameterization of the weights is different from that of
-`numpy.polyfit`, where the weights would be specified through 
-`[ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ]` 
-with the answer solving 
+`numpy.polyfit`, where the weights would be specified through
+`[ω₁,ω₂,…,ωₙ] = [√w₁, √w₂,…,√wₙ]`
+with the answer solving
 `argminᵦ | (ωᵢ⋅yᵢ- ΣⱼVᵢⱼ(ω⋅x) βⱼ) |^2`.
 
 When specified as a matrix, `W`, the solution is through the normal
@@ -131,8 +131,8 @@ fit′(P::Type{<:AbstractPolynomial},
      x::AbstractVector{T},
      y::AbstractVector{T},
      args...; kwargs...) where {T} = fit(P,x,y,args...;  kwargs...)
-         
-         
+
+
 fit(x::AbstractVector,
     y::AbstractVector,
     deg::Integer = length(x) - 1;
@@ -168,7 +168,7 @@ _wlstsq(vand, y, W::AbstractMatrix) = qr(vand' * W * vand) \ (vand' * W * y)
 """
     roots(::AbstractPolynomial; kwargs...)
 
-Returns the roots, or zeros, of the given polynomial. 
+Returns the roots, or zeros, of the given polynomial.
 
 This is calculated via the eigenvalues of the companion matrix. The `kwargs` are passed to the `LinearAlgeebra.eigvals` call.
 
@@ -181,7 +181,7 @@ This is calculated via the eigenvalues of the companion matrix. The `kwargs` are
 
 
 """
-function roots(q::AbstractPolynomial{T}; kwargs...) where {T <: Number}
+function roots(q::AbstractPolynomial{T}; kwargs...) where {T}
 
     p = convert(Polynomial{T},  q)
     roots(p; kwargs...)
@@ -231,7 +231,7 @@ end
 
 Compute the definite integral of the given polynomial from `a` to `b`. Will throw an error if either `a` or `b` are out of the polynomial's domain.
 """
-function integrate(p::AbstractPolynomial, a::Number, b::Number)
+function integrate(p::AbstractPolynomial, a, b)
     P = integrate(p)
     return P(b) - P(a)
 end
@@ -328,8 +328,8 @@ function chop!(ps::Vector{T};
     tol = norm(ps) * rtol + atol
     for i = lastindex(ps):-1:1
         val = ps[i]
-        if abs(val) > tol 
-            resize!(ps, i); 
+        if abs(val) > tol
+            resize!(ps, i);
             return nothing
         end
     end
@@ -428,8 +428,12 @@ LinearAlgebra.transpose!(p::AbstractPolynomial) = p
 Conversions =#
 Base.convert(::Type{P}, p::P) where {P <: AbstractPolynomial} = p
 Base.convert(P::Type{<:AbstractPolynomial}, x) = P(x)
-function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:Polynomials.AbstractPolynomial{T}}
-    Polynomials.isconstant(p) && return convert(S, Polynomials.constantterm(p))
+function Base.convert(::Type{S}, p::P) where {S <: Number,T, P<:AbstractPolynomial{T}}
+    isconstant(p) && return convert(S, constantterm(p))
+    throw(ArgumentError("Can't convert a nonconstant polynomial to type $S"))
+end
+function Base.convert(::Type{T}, p::P) where {T, P<:AbstractPolynomial{T}}
+    isconstant(p) && return constantterm(p)
     throw(ArgumentError("Can't convert a nonconstant polynomial to type $S"))
 end
 
@@ -537,8 +541,9 @@ Determine whether a polynomial is a monic polynomial, i.e., its leading coeffici
 ismonic(p::AbstractPolynomial) = isone(convert(Polynomial, p)[end])
 
 "`hasnan(p::AbstractPolynomial)` are any coefficients `NaN`"
-hasnan(p::AbstractPolynomial) = any(isnan, p)
-
+hasnan(p::AbstractPolynomial) = any(hasnan, p)
+hasnan(p::AbstractArray) = any(hasnan.(p))
+hasnan(x) = isnan(x)
 
 """
     isconstant(::AbstractPolynomial)
@@ -554,6 +559,8 @@ Return the coefficient vector. For a standard basis polynomial these are `[a_0,
 """
 coeffs(p::AbstractPolynomial) = p.coeffs
 
+# hook in for offset coefficients of Laurent Polynomials
+_coeffs(p::AbstractPolynomial) = coeffs(p)
 
 
 # specialize this to p[0] when basis vector is 1
@@ -570,7 +577,7 @@ constantterm(p::AbstractPolynomial{T}) where {T} = p(zero(T))
 Return the degree of the polynomial, i.e. the highest exponent in the polynomial that
 has a nonzero coefficient. The degree of the zero polynomial is defined to be -1. The default method assumes the basis polynomial, `βₖ` has degree `k`.
 """
-degree(p::AbstractPolynomial) = iszero(p) ? -1 : lastindex(p) 
+degree(p::AbstractPolynomial) = iszero(p) ? -1 : lastindex(p)
 
 
 """
@@ -617,7 +624,7 @@ Base.eachindex(p::AbstractPolynomial) = 0:length(p) - 1
 Base.broadcastable(p::AbstractPolynomial) = Ref(p)
 
 # getindex
-function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T <: Number}
+function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T}
     m,M = firstindex(p), lastindex(p)
     idx < m && throw(BoundsError(p, idx))
     idx > M && return zero(T)
@@ -628,7 +635,7 @@ Base.getindex(p::AbstractPolynomial, indices) = [p[i] for i in indices]
 Base.getindex(p::AbstractPolynomial, ::Colon) = coeffs(p)
 
 # setindex
-function Base.setindex!(p::AbstractPolynomial, value::Number, idx::Int)
+function Base.setindex!(p::AbstractPolynomial, value, idx::Int)
     n = length(coeffs(p))
     if n ≤ idx
         resize!(p.coeffs, idx + 1)
@@ -638,16 +645,18 @@ function Base.setindex!(p::AbstractPolynomial, value::Number, idx::Int)
     return p
 end
 
-Base.setindex!(p::AbstractPolynomial, value::Number, idx::Number) =
+Base.setindex!(p::AbstractPolynomial, value, idx::Number) =
     setindex!(p, value, convert(Int, idx))
-Base.setindex!(p::AbstractPolynomial, value::Number, indices) =
-    [setindex!(p, value, i) for i in indices]
+#Base.setindex!(p::AbstractPolynomial, value::Number, indices) =
+#    [setindex!(p, value, i) for i in indices]
 Base.setindex!(p::AbstractPolynomial, values, indices) =
-    [setindex!(p, v, i) for (v, i) in zip(values, indices)]
-Base.setindex!(p::AbstractPolynomial, value::Number, ::Colon) =
-    setindex!(p, value, eachindex(p))
+    [setindex!(p, v, i) for (v, i) in tuple.(values, indices)]
+#    [setindex!(p, v, i) for (v, i) in zip(values, indices)]
+#Base.setindex!(p::AbstractPolynomial, value, ::Colon) =
+#    setindex!(p, value, eachindex(p))
 Base.setindex!(p::AbstractPolynomial, values, ::Colon) =
-    [setindex!(p, v, i) for (v, i) in zip(values, eachindex(p))]
+#        [setindex!(p, v, i) for (v, i) in zip(values, eachindex(p))]
+    [setindex!(p, v, i) for (v, i) in tuple.(values, eachindex(p))]
 
 #=
 Iteration =#
@@ -826,11 +835,17 @@ basis(p::P, k::Int, _var=indeterminate(p); var=_var) where {P<:AbstractPolynomia
 
 #=
 arithmetic =#
-Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p)) 
+Base.:-(p::P) where {P <: AbstractPolynomial} = _convert(p, -coeffs(p))
+
+Base.:*(p::AbstractPolynomial, c::Number) = scalar_mult(p, c)
+Base.:*(c::Number, p::AbstractPolynomial) = scalar_mult(c, p)
+Base.:*(c::T, p::P) where {T, P <: AbstractPolynomial{T}} = scalar_mult(c, p)
+Base.:*(p::P, c::T) where {T, P <: AbstractPolynomial{T}} = scalar_mult(p, c)
+
+# implicitly identify c::Number with a constant polynomials
 Base.:+(c::Number, p::AbstractPolynomial) = +(p, c)
 Base.:-(p::AbstractPolynomial, c::Number) = +(p, -c)
 Base.:-(c::Number, p::AbstractPolynomial) = +(-p, c)
-Base.:*(c::Number, p::AbstractPolynomial) = *(p, c)
 
 # scalar operations
 # no generic p+c, as polynomial addition falls back to scalar ops
@@ -849,7 +864,7 @@ Base.:-(p1::AbstractPolynomial, p2::AbstractPolynomial) = +(p1, -p2)
 ## +(p::P,c::Number) and +(p::P, q::Q) where {T,S,X,P<:SubtypePolynomial{T,X},Q<:SubtypePolynomial{S,X}}
 ## though the default for poly+poly isn't terrible
 
-# polynomial + scalar
+# polynomial + scalar; implicit identification of c with c*one(P)
 Base.:+(p::P, c::T) where {T,X, P<:AbstractPolynomial{T,X}} = p + c * one(P)
 
 function Base.:+(p::P, c::S) where {T,X, P<:AbstractPolynomial{T,X}, S}
@@ -887,7 +902,7 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Vector{T}, p2::Vector{S}) where
     for i in 1:n2
         cs[i] += p2[i]
     end
-    
+
     return cs
 end
 
@@ -914,10 +929,10 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Dict{Int,T}, p2::Dict{Int,S}) wh
 
     R = promote_type(T,S)
     p = Dict{Int, R}()
-                     
+
 
     # this allocates in the union
-#    for i in union(eachindex(p1), eachindex(p2)) 
+#    for i in union(eachindex(p1), eachindex(p2))
 #        p[i] = p1[i] + p2[i]
 #    end
 
@@ -929,18 +944,28 @@ function ⊕(P::Type{<:AbstractPolynomial}, p1::Dict{Int,T}, p2::Dict{Int,S}) wh
             @inbounds p[i] = get(p1, i, zero(R)) + pi
         end
     end
-    
+
     return  p
 
 end
 
 ## -- multiplication
 
-## scalar *, /
-function Base.:*(p::P, c::S) where {P <: AbstractPolynomial,S}
-    _convert(p, coeffs(p) .* c)
+
+# this fall back not necessarily efficient (e.g., sparse)
+function scalar_mult(p::P, c::S) where {S, T, X, P<:AbstractPolynomial{T,X}}
+    R = Base.promote_op(*, T, S) # typeof(one(T)*one(S))?
+    𝐏 = ⟒(P){R,X}
+    𝐏([pᵢ * c for pᵢ ∈ coeffs(p)])
 end
 
+function scalar_mult(c::S, p::P) where {S, T, X, P<:AbstractPolynomial{T, X}}
+    R = Base.promote_op(*, T, S)
+    𝐏 = ⟒(P){R,X}
+    𝐏([c * pᵢ for pᵢ ∈ coeffs(p)])
+end
+
+
 function Base.:/(p::P, c::S) where {P <: AbstractPolynomial,S}
     _convert(p, coeffs(p) ./ c)
 end
@@ -955,7 +980,6 @@ function Base.:*(p1::P, p2::Q) where {T,X,P <: AbstractPolynomial{T,X},S,Y,Q <:
     return p1 * p2
 end
 
-
 Base.:^(p::AbstractPolynomial, n::Integer) = Base.power_by_squaring(p, n)
 
 function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial}
@@ -1025,7 +1049,7 @@ Comparisons =#
 Base.isequal(p1::P, p2::P) where {P <: AbstractPolynomial} = hash(p1) == hash(p2)
 Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) =
     check_same_variable(p1,p2) && (coeffs(p1) == coeffs(p2))
-Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && p[0] == n
+Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && constantterm(p) == n
 Base.:(==)(n::Number, p::AbstractPolynomial) = p == n
 
 function Base.isapprox(p1::AbstractPolynomial{T,X},
@@ -1050,7 +1074,7 @@ function Base.isapprox(p1::AbstractPolynomial{T},
                        n::S;
                        rtol::Real = (Base.rtoldefault(T, S, 0)),
                        atol::Real = 0,) where {T,S}
-    return isapprox(p1, _convert(p1, [n])) 
+    return isapprox(p1, _convert(p1, [n]))
 end
 
 Base.isapprox(n::S,
diff --git a/src/contrib.jl b/src/contrib.jl
index a5d20b7c..9897ff3a 100644
--- a/src/contrib.jl
+++ b/src/contrib.jl
@@ -2,6 +2,7 @@
 
 using Base.Cartesian
 
+
 # direct version (do not check if threshold is satisfied)
 @generated function fastconv(E::Array{T,N}, k::Array{T,N}) where {T,N}
     quote
@@ -40,7 +41,7 @@ function evalpoly(x::S, p::Tuple) where {S}
         N = length(p.parameters)
         ex = :(p[end]*_one(S))
         for i in N-1:-1:1
-            ex = :(_muladd(x, $ex, p[$i]))
+            ex = :(_muladd($ex, x, p[$i]))
         end
         ex
     else
@@ -54,7 +55,7 @@ function _evalpoly(x::S, p) where {S}
     N = length(p)
     ex = p[end]*_one(x)
     for i in N-1:-1:1
-        ex = _muladd(x, ex, p[i])
+        ex = _muladd(ex, x, p[i])
     end
     ex
 end
@@ -110,8 +111,8 @@ end
 
 ## modify muladd, as needed
 _muladd(a,b,c) = muladd(a,b,c)
-_muladd(a::Vector, b, c) = a.*b .+ c
-_muladd(a::Matrix, b, c) = a*(b*I) + c*I
+_muladd(a, b::Vector, c) = a.*b .+ c
+_muladd(a, b::Matrix, c) = (a*I)*b + c*I
 
 # try to get y = P(c::T)(x::S) = P{T}(c)(x::S) to
 # have y == one(T)*one(S)*x
@@ -129,4 +130,3 @@ end
 @generated function constructorof(::Type{T}) where T
     getfield(parentmodule(T), nameof(T))
 end
-
diff --git a/src/polynomials/ChebyshevT.jl b/src/polynomials/ChebyshevT.jl
index 1ab386c9..b5a4e511 100644
--- a/src/polynomials/ChebyshevT.jl
+++ b/src/polynomials/ChebyshevT.jl
@@ -1,7 +1,7 @@
 export ChebyshevT
 
 """
-    ChebyshevT{T<:Number, X}(coeffs::AbstractVector)
+    ChebyshevT{T, X}(coeffs::AbstractVector)
 
 Chebyshev polynomial of the first kind.
 
@@ -9,7 +9,7 @@ Construct a polynomial from its coefficients `coeffs`, lowest order first, optio
 terms of the given variable `var`, which can be a character, symbol, or string.
 
 !!! note
-    `ChebyshevT` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first 
+    `ChebyshevT` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
     index always corresponding to the coefficient of `T_0(x)`.
 
 # Examples
@@ -39,9 +39,9 @@ The latter shows how to evaluate a `ChebyshevT` polynomial outside of its domain
     The Chebyshev polynomials are also implemented in `ApproxFun`, `ClassicalOrthogonalPolynomials.jl`, `FastTransforms.jl`, and `SpecialPolynomials.jl`.
 
 """
-struct ChebyshevT{T <: Number, X} <: AbstractPolynomial{T, X}
+struct ChebyshevT{T, X} <: AbstractPolynomial{T, X}
     coeffs::Vector{T}
-    function ChebyshevT{T, X}(coeffs::AbstractVector{S}) where {T <: Number,X, S}
+    function ChebyshevT{T, X}(coeffs::AbstractVector{S}) where {T, X, S}
 
         if Base.has_offset_axes(coeffs)
             @warn "ignoring the axis offset of the coefficient vector"
@@ -264,8 +264,12 @@ end
 function showterm(io::IO, ::Type{ChebyshevT{T,X}}, pj::T, var, j, first::Bool, mimetype) where {N, T,X}
     iszero(pj) && return false
     !first &&  print(io, " ")
-    print(io, hasneg(T)  && isneg(pj) ? "- " :  (!first ? "+ " : ""))
-    print(io, "$(abs(pj))⋅T_$j($var)")
+    if hasneg(T)
+        print(io, isneg(pj) ? "- " :  (!first ? "+ " : ""))
+        print(io, "$(abs(pj))⋅T_$j($var)")
+    else
+        print(io, "+ ", "$(pj)⋅T_$j($var)")
+    end
     return true
 end
 
diff --git a/src/polynomials/ImmutablePolynomial.jl b/src/polynomials/ImmutablePolynomial.jl
index 5e81deec..1acfe344 100644
--- a/src/polynomials/ImmutablePolynomial.jl
+++ b/src/polynomials/ImmutablePolynomial.jl
@@ -1,9 +1,9 @@
 export ImmutablePolynomial
 
 """
-    ImmutablePolynomial{T<:Number, X, N}(coeffs::AbstractVector{T})
+    ImmutablePolynomial{T, X, N}(coeffs::AbstractVector{T})
 
-Construct an immutable (static) polynomial from its coefficients 
+Construct an immutable (static) polynomial from its coefficients
 `a₀, a₁, …, aₙ`,
 lowest order first, optionally in terms of the given variable `x`
 where `x` can be a character, symbol, or string.
@@ -29,7 +29,7 @@ immutable polynomials can not promote to a common type. As such, they
 are precluded from use in rational functions.
 
 !!! note
-    `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first 
+    `ImmutablePolynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
     index always corresponding to the constant term.
 
 # Examples
@@ -51,9 +51,9 @@ ImmutablePolynomial(1.0)
     This was modeled after [StaticUnivariatePolynomials](https://github.com/tkoolen/StaticUnivariatePolynomials.jl) by `@tkoolen`.
 
 """
-struct ImmutablePolynomial{T <: Number, X, N} <: StandardBasisPolynomial{T, X}
+struct ImmutablePolynomial{T, X, N} <: StandardBasisPolynomial{T, X}
     coeffs::NTuple{N, T}
-    function ImmutablePolynomial{T,X,N}(coeffs::NTuple{N,T}) where {T <: Number, X, N}
+    function ImmutablePolynomial{T,X,N}(coeffs::NTuple{N,T}) where {T, X, N}
         N == 0 && return new{T,X, 0}(coeffs)
         iszero(coeffs[end]) &&  throw(ArgumentError("Leading  term must  be  non-zero"))
         new{T,X,N}(coeffs)
@@ -63,19 +63,29 @@ end
 @register ImmutablePolynomial
 
 ## Various interfaces
-
 ## Abstract Vector coefficients
+function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{T})  where {T,X, N}
+    cs = NTuple{N,T}(coeffs[i] for i ∈ firstindex(coeffs):N)
+    ImmutablePolynomial{T, X, N}(cs)
+
+end
+
+function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{S})  where {T,X, N, S}
+    cs = NTuple{N,T}(coeffs[i] for i ∈ firstindex(coeffs):N)
+    ImmutablePolynomial{T, X, N}(cs)
+
+end
+
 function ImmutablePolynomial{T,X}(coeffs::AbstractVector{S})  where {T,X,S}
     R = promote_type(T,S)
-    
+
     if Base.has_offset_axes(coeffs)
         @warn "ignoring the axis offset of the coefficient vector"
     end
     N = findlast(!iszero, coeffs)
     N == nothing && return ImmutablePolynomial{R,X,0}(())
     N′ = N + 1 - firstindex(coeffs)
-    cs = NTuple{N′,T}(coeffs[i] for i ∈ firstindex(coeffs):N)
-    ImmutablePolynomial{T, X, N′}(cs)
+    ImmutablePolynomial{T, X, N′}([coeffs[i] for i ∈ firstindex(coeffs):N])
 end
 
 ## -- Tuple arguments
@@ -105,7 +115,7 @@ Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p))
 degree(p::ImmutablePolynomial{T,X, N}) where {T,X,N} = N - 1 # no trailing zeros
 isconstant(p::ImmutablePolynomial{T,X,N}) where {T,X,N}  = N <= 1
 
-Base.setindex!(p::ImmutablePolynomial, val::Number,  idx::Int) = throw(ArgumentError("ImmutablePolynomials are immutable"))
+Base.setindex!(p::ImmutablePolynomial, val,  idx::Int) = throw(ArgumentError("ImmutablePolynomials are immutable"))
 
 for op in [:isequal, :(==)]
     @eval function Base.$op(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
@@ -116,7 +126,7 @@ for op in [:isequal, :(==)]
         n2 = findlast(!iszero, p2s)
         (n1 == nothing && n2 == nothing) && return true
         (n1 == nothing || n2  == nothing) && return false
-        $op(p1s[1:n1],p2s[1:n2]) &&  return true 
+        $op(p1s[1:n1],p2s[1:n2]) &&  return true
         false
     end
 end
@@ -154,7 +164,7 @@ function Base.:+(p::P, c::S) where {T, X, N, P <: ImmutablePolynomial{T,X,N}, S<
     q = ImmutablePolynomial{R,X,N}(cs)
 
     return q
-    
+
 end
 
 Base.:-(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = ImmutablePolynomial{T,X,N}(.-p.coeffs)
@@ -179,15 +189,21 @@ end
 
 ## multiplication
 
-function Base.:*(p::ImmutablePolynomial{T,X,N}, c::S) where {T, X,N, S <: Number}
-    R = eltype(one(T)*one(S))
-    P = ImmutablePolynomial{R,X}
-    (N == 0  || iszero(c)) && return zero(P)
+function scalar_mult(p::ImmutablePolynomial{T,X,N}, c::S) where {T, X,N, S <: Number}
+    R = eltype(p[0] * c * 0)
+    (N == 0  || iszero(c)) && return zero(ImmutablePolynomial{R,X})
     cs = p.coeffs .* c
-    iszero(cs[end]) ? P(cs) : P{N}(cs) # more performant to specify when N is known
+    return ImmutablePolynomial(cs, X)
+ end
+
+function scalar_mult(c::S, p::ImmutablePolynomial{T,X,N}) where {T, X,N, S <: Number}
+    R = eltype(p[0] * c * 0)
+    (N == 0  || iszero(c)) && return zero(ImmutablePolynomial{R,X})
+    cs = p.coeffs .* c
+    return ImmutablePolynomial(cs, X)
 end
 
-function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S <: Number}
+function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S}
     R = eltype(one(T)/one(S))
     P = ImmutablePolynomial{R,X}
     (N == 0  || isinf(c)) && return zero(P)
@@ -222,12 +238,12 @@ LinearAlgebra.norm(q::ImmutablePolynomial, p::Real) = _norm(q.coeffs, p)
     for j = 2:N
         expr = :($expr + abs2(a[$j]))
     end
-    
+
     return quote
         $(Expr(:meta, :inline))
         @inbounds return sqrt($expr)
     end
-    
+
 end
 
 _norm_p0(x) = iszero(x) ? zero(x) : one(x)
@@ -258,5 +274,3 @@ _norm_p0(x) = iszero(x) ? zero(x) : one(x)
     end
 
 end
-
-
diff --git a/src/polynomials/LaurentPolynomial.jl b/src/polynomials/LaurentPolynomial.jl
index 31403516..4527dc6a 100644
--- a/src/polynomials/LaurentPolynomial.jl
+++ b/src/polynomials/LaurentPolynomial.jl
@@ -5,7 +5,7 @@ export LaurentPolynomial
 
 A [Laurent](https://en.wikipedia.org/wiki/Laurent_polynomial) polynomial is of the form `a_{m}x^m + ... + a_{n}x^n` where `m,n` are  integers (not necessarily positive) with ` m <= n`.
 
-The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`. 
+The `coeffs` specify `a_{m}, a_{m-1}, ..., a_{n}`.
 The argument `m` represents the lowest exponent of the variable in the series, and is taken to be zero by default.
 
 Laurent polynomials and standard basis polynomials promote to  Laurent polynomials. Laurent polynomials may be  converted to a standard basis  polynomial when `m >= 0`
@@ -69,12 +69,12 @@ julia> x^degree(p) * p(x⁻¹) # reverses  coefficients
 LaurentPolynomial(3.0 + 2.0*x + 1.0*x²)
 ```
 """
-struct LaurentPolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
+struct LaurentPolynomial{T, X} <: StandardBasisPolynomial{T, X}
     coeffs::Vector{T}
     m::Base.RefValue{Int}
     n::Base.RefValue{Int}
     function LaurentPolynomial{T,X}(coeffs::AbstractVector{S},
-                                    m::Union{Int, Nothing}=nothing) where {T <: Number, X, S}
+                                    m::Union{Int, Nothing}=nothing) where {T, X, S}
 
         fnz = findfirst(!iszero, coeffs)
         fnz == nothing && return  new{T,X}(zeros(T,1), Ref(0), Ref(0))
@@ -84,7 +84,7 @@ struct LaurentPolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
             cs = convert(Vector{T}, coeffs[fnz:lnz])
             return new{T,X}(cs, Ref(fnz), Ref(lnz))
         else
-            
+
             c = convert(Vector{T}, coeffs[fnz:lnz])
 
             m′ = fnz - 1 + (m == nothing ? 0 : m)
@@ -100,16 +100,15 @@ end
 
 ## constructors
 function LaurentPolynomial{T}(coeffs::AbstractVector{S}, m::Int, var::SymbolLike=:x) where {
-    T <: Number, S <: Number}
+    T, S <: Number}
     LaurentPolynomial{T,Symbol(var)}(T.(coeffs), m)
 end
 
-function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {
-    T <: Number}
+function LaurentPolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T}
     LaurentPolynomial{T, Symbol(var)}(coeffs, 0)
 end
 
-function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T <: Number}
+function LaurentPolynomial(coeffs::AbstractVector{T}, m::Int, var::SymbolLike=:x) where {T}
     LaurentPolynomial{T, Symbol(var)}(coeffs, m)
 end
 
@@ -133,21 +132,17 @@ function Base.convert(P::Type{<:Polynomial}, q::LaurentPolynomial)
     P([q[i] for i  in 0:n], indeterminate(q))
 end
 
-# need to add p.m[], so abstract.jl method isn't sufficent
-function Base.convert(::Type{P}, p::LaurentPolynomial) where {P<:LaurentPolynomial}
-    S′ = _eltype(P)
-    Y′ = _indeterminate(P)
-    S = S′ == nothing ? eltype(p) : S′
-    Y = Y′ == nothing ? indeterminate(p) : Y′
-    isconstant(p) && return LaurentPolynomial{S,Y}(constantterm(p))
-    LaurentPolynomial{S,Y}(p.coeffs, p.m[])
+
+# work around for non-applicable convert(::Type{<:P}, p::P{T,X}) in abstract.jl
+struct OffsetCoeffs{V}
+    coeffs::V
+    m::Int
 end
 
-# function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {T, P <:LaurentPolynomial{T},S}
-#     v′ = _indeterminate(P)
-#     X = v′ == nothing ? indeterminate(q) : v′
-#     ⟒(P){T,X}([q[i] for i in 0:degree(q)], 0)
-# end
+_coeffs(p::LaurentPolynomial) = OffsetCoeffs(p.coeffs, p.m[])
+function LaurentPolynomial{T,X}(p::OffsetCoeffs) where {T, X}
+    LaurentPolynomial{T,X}(p.coeffs, p.m)
+end
 
 function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {P <:LaurentPolynomial,S}
 
@@ -528,7 +523,7 @@ function integrate(p::P) where {T, X, P<: LaurentPolynomial{T, X}}
     !iszero(p[-1])  && throw(ArgumentError("Can't integrate Laurent  polynomial with  `x⁻¹` term"))
     R = eltype(one(T)/1)
     Q = ⟒(P){R, X}
-    
+
     if hasnan(p)
         return Q([NaN],0)
     end
diff --git a/src/polynomials/Polynomial.jl b/src/polynomials/Polynomial.jl
index 4a5ca6b8..5cb46223 100644
--- a/src/polynomials/Polynomial.jl
+++ b/src/polynomials/Polynomial.jl
@@ -1,7 +1,7 @@
 export Polynomial
 
 """
-    Polynomial{T<:Number, X}(coeffs::AbstractVector{T}, [var = :x])
+    Polynomial{T, X}(coeffs::AbstractVector{T}, [var = :x])
 
 Construct a polynomial from its coefficients `coeffs`, lowest order first, optionally in
 terms of the given variable `var` which may be a character, symbol, or a string.
@@ -14,8 +14,8 @@ with combinations of polynomials and scalars. However, operations involving two
 polynomials of different variables causes an error except those involving a constant polynomial.
 
 !!! note
-    `Polynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first 
-    index always corresponding to the constant term. In order to use the axis of `coeffs` as exponents, 
+    `Polynomial` is not axis-aware, and it treats `coeffs` simply as a list of coefficients with the first
+    index always corresponding to the constant term. In order to use the axis of `coeffs` as exponents,
     consider using a [`LaurentPolynomial`](@ref) or possibly a [`SparsePolynomial`](@ref).
 
 # Examples
@@ -32,9 +32,9 @@ julia> one(Polynomial)
 Polynomial(1.0)
 ```
 """
-struct Polynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
+struct Polynomial{T, X} <: StandardBasisPolynomial{T, X}
     coeffs::Vector{T}
-    function Polynomial{T, X}(coeffs::AbstractVector{S}) where {T <: Number, X, S}
+    function Polynomial{T, X}(coeffs::AbstractVector{S}) where {T, X, S}
         if Base.has_offset_axes(coeffs)
             @warn "ignoring the axis offset of the coefficient vector"
         end
@@ -44,7 +44,7 @@ struct Polynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
         new{T,X}(cs)
     end
     # non-copying alternative assuming !iszero(coeffs[end])
-    function Polynomial{T, X}(checked::Val{false}, coeffs::AbstractVector{T}) where {T <: Number, X}
+    function Polynomial{T, X}(checked::Val{false}, coeffs::AbstractVector{T}) where {T, X}
         new{T, X}(coeffs)
     end
 end
@@ -56,7 +56,7 @@ end
 
 # scalar +,* faster  than standard-basis/common versions as it avoids a copy
 function Base.:+(p::P, c::S) where {T, X, P <: Polynomial{T, X}, S<:Number}
-    R = promote_type(T, S)
+    R = Base.promote_op(+, T, S)
     Q = Polynomial{R,X}
     as = convert(Vector{R}, copy(coeffs(p)))
     as[1] += c
@@ -64,7 +64,7 @@ function Base.:+(p::P, c::S) where {T, X, P <: Polynomial{T, X}, S<:Number}
 end
 
 function Base.:*(p::P, c::S) where {T, X, P <: Polynomial{T,X} , S <: Number}
-    R = promote_type(T,S)
+    R = Base.promote_op(*, T, S) #promote_type(T,S)
     Q = Polynomial{R, X}
     as = R[aᵢ * c for aᵢ ∈ coeffs(p)]
     iszero(as[end]) ? Q(as) : Q(Val(false), as)
@@ -88,9 +88,9 @@ function Base.:+(p1::P1, p2::P2) where {T,X, P1<:Polynomial{T,X},
 
     Q(Val(false), cs)
 end
-        
-function Base.:*(p::P, q::P) where {T,X, P<:Polynomial{T,X}}
+
+# redundant, a bit faster
+function Base.:*(p::P, q::P) where {T <: Number,X, P<:Polynomial{T,X}}
     c = fastconv(p.coeffs, q.coeffs)
     return iszero(c[end]) ? P(c) : P(Val(false), c)
 end
-
diff --git a/src/polynomials/SparsePolynomial.jl b/src/polynomials/SparsePolynomial.jl
index a2e0a69f..aa7995e9 100644
--- a/src/polynomials/SparsePolynomial.jl
+++ b/src/polynomials/SparsePolynomial.jl
@@ -5,7 +5,7 @@ export   SparsePolynomial
 
 Polynomials in the standard basis backed by a dictionary holding the
 non-zero coefficients. For polynomials of high degree, this might be
-advantageous. 
+advantageous.
 
 # Examples:
 
@@ -38,31 +38,31 @@ julia> p(1)
 ```
 
 """
-struct SparsePolynomial{T <: Number, X} <: StandardBasisPolynomial{T, X}
+struct SparsePolynomial{T, X} <: StandardBasisPolynomial{T, X}
     coeffs::Dict{Int, T}
-    function SparsePolynomial{T, X}(coeffs::AbstractDict{Int, S}) where {T <: Number, X, S}
+    function SparsePolynomial{T, X}(coeffs::AbstractDict{Int, S}) where {T, X, S}
         c = Dict{Int, T}(coeffs)
         for (k,v)  in coeffs
             iszero(v) && pop!(c,  k)
         end
         new{T, X}(c)
     end
-    function SparsePolynomial{T,X}(checked::Val{false}, coeffs::AbstractDict{Int, T}) where {T <: Number, X}
+    function SparsePolynomial{T,X}(checked::Val{false}, coeffs::AbstractDict{Int, T}) where {T, X}
         new{T,X}(coeffs)
     end
 end
 
 @register SparsePolynomial
 
-function SparsePolynomial{T}(coeffs::AbstractDict{Int, S}, var::SymbolLike=:x) where {T <: Number, S}
+function SparsePolynomial{T}(coeffs::AbstractDict{Int, S}, var::SymbolLike=:x) where {T, S}
     SparsePolynomial{T, Symbol(var)}(convert(Dict{Int,T}, coeffs))
 end
 
-function SparsePolynomial(coeffs::AbstractDict{Int, T}, var::SymbolLike=:x) where {T <: Number}
+function SparsePolynomial(coeffs::AbstractDict{Int, T}, var::SymbolLike=:x) where {T}
     SparsePolynomial{T, Symbol(var)}(coeffs)
 end
 
-function SparsePolynomial{T,X}(coeffs::AbstractVector{S}) where {T <: Number, X, S}
+function SparsePolynomial{T,X}(coeffs::AbstractVector{S}) where {T, X, S}
 
     if Base.has_offset_axes(coeffs)
         @warn "ignoring the axis offset of the coefficient vector"
@@ -109,11 +109,11 @@ function coeffs(p::SparsePolynomial{T})  where {T}
         cs[k+1]=v
     end
     cs
-    
+
 end
 
 # get/set index
-function Base.getindex(p::SparsePolynomial{T}, idx::Int) where {T <: Number}
+function Base.getindex(p::SparsePolynomial{T}, idx::Int) where {T}
     get(p.coeffs, idx, zero(T))
 end
 
@@ -151,16 +151,16 @@ end
 ##
 ## ----
 ##
-    
+
 function evalpoly(x::S, p::SparsePolynomial{T}) where {T,S}
 
-    tot = zero(T) * EvalPoly._one(x) 
+    tot = zero(T) * EvalPoly._one(x)
     for (k,v) in p.coeffs
         tot = EvalPoly._muladd(x^k, v, tot)
     end
-    
+
     return tot
-    
+
 end
 
 #  map: over values -- not keys
@@ -186,7 +186,7 @@ function Base.:+(p::SparsePolynomial{T,X}, c::S) where {T, X, S <: Number}
     iszero(D[0]) && pop!(D,0)
 
     return P(Val(false), D)
-    
+
 end
 
 # Implement over fallback. A bit faster as it covers T != S
@@ -195,28 +195,42 @@ function Base.:+(p1::P1, p2::P2) where {T,X, P1<:SparsePolynomial{T,X},
 
     R = promote_type(T,S)
     Q = SparsePolynomial{R,X}
-    
+
     d1, d2 = degree(p1), degree(p2)
     cs = d1 > d2 ? ⊕(P1, p1.coeffs, p2.coeffs) : ⊕(P1, p2.coeffs, p1.coeffs)
 
     return d1 != d2 ? Q(Val(false), cs) : Q(cs)
-    
+
 end
 
 ## Multiplication
-function Base.:*(p::P, c::S) where {T, X, P <: SparsePolynomial{T,X}, S <: Number}
+function scalar_mult(p::P, c::S) where {T, X, P <: SparsePolynomial{T,X}, S<:Number}
 
-    R = promote_type(T,S)
+    R1 = promote_type(T,S)
+    R = typeof(zero(c)*zero(T))
     Q = ⟒(P){R,X}
-    
-    q  = zero(Q)
+    q = zero(Q)
     for (k,pₖ) ∈ pairs(p)
         q[k] = pₖ * c
     end
-    
+
     return q
 end
 
+function scalar_mult(c::S, p::P) where {T, X, P <: SparsePolynomial{T,X}, S<:Number}
+
+    R1 = promote_type(T,S)
+    R = typeof(zero(c)*zero(T))
+    Q = ⟒(P){R,X}
+    q = zero(Q)
+    for (k,pₖ) ∈ pairs(p)
+        q[k] = c * pₖ
+    end
+
+    return q
+end
+
+
 function Base.:*(p::P, q::Q) where {T,X,P<:SparsePolynomial{T,X},
                                     S,  Q<:SparsePolynomial{S,X}}
     R = promote_type(T,S)
@@ -225,7 +239,7 @@ end
 
 
 function derivative(p::SparsePolynomial{T,X}, order::Integer = 1) where {T,X}
-    
+
     order < 0 && throw(ArgumentError("Order of derivative must be non-negative"))
     order == 0 && return p
 
@@ -247,7 +261,7 @@ end
 
 
 function integrate(p::P) where {T, X, P<:SparsePolynomial{T,X}}
-    
+
     R = eltype(one(T)/1)
     Q = SparsePolynomial{R,X}
 
@@ -263,4 +277,3 @@ function integrate(p::P) where {T, X, P<:SparsePolynomial{T,X}}
     return ∫p
 
 end
-
diff --git a/src/polynomials/factored_polynomial.jl b/src/polynomials/factored_polynomial.jl
index 25dbe067..3f4c07db 100644
--- a/src/polynomials/factored_polynomial.jl
+++ b/src/polynomials/factored_polynomial.jl
@@ -142,7 +142,7 @@ function printpoly(io::IO, p::FactoredPolynomial{T,X}, mimetype=nothing) where {
                 print(io, "(")
                 print(io, x)
                 if hasneg(T)
-                    isneg(k) ?  print(io, " + ", -k) : print(io, " - ", k) 
+                    isneg(k) ?  print(io, " + ", -k) : print(io, " - ", k)
                 else
                     print(io, " - ", k)
                 end
@@ -191,13 +191,13 @@ function Base.isapprox(p1::FactoredPolynomial{T,X},
     # # sorting roots below works only with real roots...
     # isapprox(p1.c, p2.c, rtol=rtol, atol=atol) || return false
     # k1,k2 = sort(collect(keys(p1.coeffs)),by = x -> (real(x), imag(x))), sort(collect(keys(p2.coeffs)),by = x -> (real(x), imag(x)))
-    
+
     # length(k1) == length(k2) || return false
     # for (k₁, k₂) ∈ zip(k1, k2)
     #     isapprox(k₁, k₂, atol=atol, rtol=rtol) || return false
     #     p1.coeffs[k₁] == p2.coeffs[k₂] || return false
     # end
-    
+
     # return true
 end
 
@@ -247,7 +247,7 @@ roots(p::FactoredPolynomial{T}) where {T} = Base.typed_vcat(T,[repeat([k],v) for
 # unary subtraction
 Base.:-(p::P) where {T,X,P<:FactoredPolynomial{T,X}} = (-1)*p
 
-# addition 
+# addition
 function Base.:+(p::P, q::P) where {T,X,P<:FactoredPolynomial{T,X}}
     𝑷 = Polynomial{T,X}
     𝒑,𝒒 = convert(𝑷, p), convert(𝑷, q)
@@ -264,7 +264,7 @@ function Base.:*(p::P, q::P) where {T,X, P<:FactoredPolynomial{T,X}}
 end
 
 # scalar mult
-function Base.:*(p::P, c::S) where {S<:Number, T, X, P <: FactoredPolynomial{T, X}}
+function scalar_mult(p::P, c::S) where {S<:Number, T, X, P <: FactoredPolynomial{T, X}}
     R = promote_type(T,S)
     d = Dict{R, Int}() # wident
     copy!(d, p.coeffs)
@@ -304,7 +304,7 @@ function uvw(p::P, q::P; kwargs...) where {T, X, P<:FactoredPolynomial{T,X}}
     du, dv, dw = Dict{T,Int}(), Dict{T,Int}(), Dict{T,Int}()
     dp,dq = p.coeffs, q.coeffs
     kp,kq = keys(dp), keys(dq)
-    
+
     for k ∈ setdiff(kp, kq)
         dv[k] = dp[k]
     end
@@ -315,7 +315,7 @@ function uvw(p::P, q::P; kwargs...) where {T, X, P<:FactoredPolynomial{T,X}}
         pₖ,qₖ = dp[k], dq[k]
         m = min(pₖ, qₖ)
         du[k] = m
-        dv[k] = pₖ - m; 
+        dv[k] = pₖ - m;
         dw[k] = qₖ - m
     end
     P(du), P(dv, p.c), P(dw, q.c)
@@ -332,7 +332,7 @@ function Base.divrem(p::P, q::P) where {T, X, P<:FactoredPolynomial{T,X}}
 
     return (d, u*r)
 
-end    
+end
 
 ## ----
 
@@ -349,4 +349,3 @@ function derivative(p::P,n::Int) where {P <: FactoredPolynomial}
     𝒑⁽ⁿ⁾ = derivative(𝒑, n)
     convert(P, 𝒑⁽ⁿ⁾)
 end
-
diff --git a/src/polynomials/standard-basis.jl b/src/polynomials/standard-basis.jl
index 0023b91b..20e08e1c 100644
--- a/src/polynomials/standard-basis.jl
+++ b/src/polynomials/standard-basis.jl
@@ -64,6 +64,9 @@ function variable(::Type{P}) where {P<:StandardBasisPolynomial}
     ⟒(P){T,X}([zero(T),one(T)])
 end
 
+# can bypass c*one(P)
+Base.:+(p::P, c::T) where {T, P<:StandardBasisPolynomial{T}} = p + ⟒(P)([c])
+
 ## multiplication algorithms for computing p * q.
 ## default multiplication between same type.
 ## subtypes might relax to match T,S to avoid one conversion
@@ -71,13 +74,32 @@ function Base.:*(p::P, q::P) where {T,X, P<:StandardBasisPolynomial{T,X}}
     cs = ⊗(P, coeffs(p), coeffs(q))
     P(cs)
 end
-                                    
-## put here, not with type defintion, in case reuse is possible
-function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T,S}
+
+function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T<:Number,S<:Number}
     R = promote_type(T,S)
     fastconv(convert(Vector{R}, p), convert(Vector{R},q))
 end
 
+## put here, not with type defintion, in case reuse is possible
+## `conv` can be used with matrix entries, unlike `fastconv`
+function conv(p::Vector{T}, q::Vector{S}) where {T,S}
+        as = [p[1]*q[1]]
+    z = 0 * as[1]
+    n,m = length(p)-1, length(q)-1
+    for i ∈ 1:n+m
+        Σ = z
+        for j ∈ max(0, i-m):min(i,n)
+            Σ += p[1+j]*q[1 + i-j]
+        end
+        push!(as, Σ)
+    end
+    as
+end
+
+function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T,S}
+    conv(p, q)
+end
+
 ## Static size of product makes generated functions  a good choice
 ## from https://github.com/tkoolen/StaticUnivariatePolynomials.jl/blob/master/src/monomial_basis.jl
 ## convolution of two tuples
@@ -97,7 +119,7 @@ end
 
     return quote
         Base.@_inline_meta
-        tuple($(exprs...))        
+        tuple($(exprs...))
     end
 
 end
@@ -129,6 +151,7 @@ function fromroots(P::Type{<:StandardBasisPolynomial}, r::AbstractVector{T}; var
     return P(reverse(c), var)
 end
 
+## ----
 
 function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPolynomial{T, X}}
     order < 0 && throw(ArgumentError("Order of derivative must be non-negative"))
@@ -136,9 +159,9 @@ function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPol
     # we avoid usage like Base.promote_op(*, T, Int) here, say, as
     # Base.promote_op(*, Rational, Int) is Any, not Rational in analogy to
     # Base.promote_op(*, Complex, Int)
-    R = eltype(one(T)*1)
+    R = typeof(constantterm(p)*1)
     Q = ⟒(P){R,X}
-    
+
     order == 0 && return p
     hasnan(p) && return Q(R[NaN])
     order > length(p) && return zero(Q)
@@ -153,15 +176,16 @@ function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPol
 end
 
 function integrate(p::P) where {T, X, P <: StandardBasisPolynomial{T, X}}
-    R = eltype(one(T)/1)
-    Q = ⟒(P){R,X}    
+
+    R = typeof(constantterm(p)/1)
+    Q = ⟒(P){R,X}
 
     hasnan(p) && return Q([NaN])
     iszero(p) && return zero(Q)
 
     n = length(p)
     as = Vector{R}(undef, n + 1)
-    as[1] = zero(R)
+    as[1] = 0*constantterm(p)
     for (i, pᵢ) ∈ pairs(p)
         i′ = i + 1
         @inbounds as[i′+1] = pᵢ/i′
@@ -171,7 +195,7 @@ end
 
 function Base.divrem(num::P, den::Q) where {T, P <: StandardBasisPolynomial{T}, S, Q <: StandardBasisPolynomial{S}}
 
-    assert_same_variable(num, den) 
+    assert_same_variable(num, den)
     X = indeterminate(num)
 
 
@@ -363,6 +387,52 @@ function companion(p::P) where {T, P <: StandardBasisPolynomial{T}}
     return comp
 end
 
+# block companion matrix
+function companion(p::P) where {T, M <: AbstractMatrix{T}, P<:StandardBasisPolynomial{M}}
+    C₀, C₁ = companion_pencil(p)
+    C₀ * inv(C₁) # could be more efficient
+end
+
+# the companion pencil is `C₀`, `C₁` where `λC₁ - C₀` is singular for
+# eigen values of the companion matrix: `vᵀ(λC₁ - C₀) = 0` => `vᵀλ = vᵀ(C₀C₁⁻¹)`, where `C₀C₁⁻¹`
+# is the companion matrix.
+function companion_pencil(p::P) where {T, P<:StandardBasisPolynomial{T}}
+    n = degree(p)
+    C₁ = diagm(0 => ones(T, n))
+    C₁[end,end] = p[end]
+
+    C₀ = diagm(-1 => ones(T, n-1))
+    for i ∈ 0:n-1
+        C₀[1+i,end] = -p[i]
+    end
+    C₀, C₁
+end
+
+# block version
+function companion_pencil(p::P) where {T, M <: AbstractMatrix{T}, P<:StandardBasisPolynomial{M}}
+
+    m,m′ = size(p[0])
+    @assert m == m′  # square matrix
+
+    n = degree(p)
+
+    C₀ = zeros(T, n*m, n*m)
+    C₁ = zeros(T, n*m, n*m)
+
+
+    Δ = 1:m
+    for i ∈ 1:n-1
+        C₁[(i-1)*m .+ Δ, (i-1)*m .+ Δ] .= I(m)
+        C₀[i*m .+ Δ, (i-1)*m .+ Δ] .= I(m)
+    end
+    C₁[(n-1)*m .+ Δ, (n-1)*m .+ Δ] .= p[end]
+    for i ∈ 0:n-1
+        C₀[i*m .+ Δ, (n-1)*m .+ Δ] .= -p[i]
+    end
+
+    C₀, C₁
+end
+
 function  roots(p::P; kwargs...)  where  {T, P <: StandardBasisPolynomial{T}}
 
     R = eltype(one(T)/one(T))
@@ -381,7 +451,7 @@ function  roots(p::P; kwargs...)  where  {T, P <: StandardBasisPolynomial{T}}
 
     k  == K && return zeros(R, k-1)
 
-    # find eigenvalues of the  companion matrix
+    # find eigenvalues of the companion matrix of the 0-deflated polynomial
     comp  = companion(⟒(P)(as[k:K], indeterminate(p)))
     L = eigvals(comp; kwargs...)
     append!(L, zeros(eltype(L), k-1))
diff --git a/src/show.jl b/src/show.jl
index 1992472b..7c716e0f 100644
--- a/src/show.jl
+++ b/src/show.jl
@@ -7,7 +7,10 @@ export printpoly
 ## which can be modified by users for other Ts
 
 _iszero(x::T) where {T} = (x == zero(T)) === true
+_iszero(x::AbstractArray{T}) where {T} = all(isequal.(x, zero(T)))
+
 _isone(x::T) where {T} = (x == one(T)) === true
+_isone(x::AbstractArray{T}) where {T} = all(isequal.(x, one(T)))
 
 "`hasneg(::T)` attribute is true if: `pj < zero(T)` is defined."
 hasneg(::Type{T}) where {T} = false
@@ -216,7 +219,7 @@ julia> Polynomial([Dual(1,2), Dual(3,4)])
 Polynomial(1 + 2ɛ + 3 + 4ɛ*x)
 ```
 
-```jldoctest 
+```jldoctest
 julia> using DualNumbers, Polynomials
 
 julia> function Base.show_unquoted(io::IO, pj::Dual, indent::Int, prec::Int)
diff --git a/test/Project.toml b/test/Project.toml
index abc0a286..0f1df4db 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,4 +5,5 @@ OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/StandardBasis.jl b/test/StandardBasis.jl
index 1c8c1778..a83b6c5f 100644
--- a/test/StandardBasis.jl
+++ b/test/StandardBasis.jl
@@ -60,7 +60,6 @@ isimmutable(::Type{<:ImmutablePolynomial}) = true
     end
 
 end
-        
 
 @testset "Mapdomain" begin
     for P in Ps
@@ -141,11 +140,220 @@ Base.getindex(z::ZVector, I::Int) = parent(z)[I + z.offset]
     end
 end
 
+if VERSION >= v"1.5"
+@testset "Non-number type" begin
+    conv = Polynomials.conv
+    @testset "T=Polynomial{Int,:y}" begin
+        for P in (Polynomial,)
+
+            T = P{Int, :y}
+            a,b,c = T([1]), T([1,2]), T([1,2,3])
+            p = P([a,b,c])
+            q = P([a,b])
+            s = 2
+            d = c
+
+            # scalar product
+            @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]])
+            @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]])
+            @test_throws ArgumentError d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]]) # can't fix
+            @test_throws ArgumentError p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]]) # can't fix
+
+            # poly add
+            @test p + q == P([a+a,b+b,c])
+            @test p - q == P([a-a,b-b,c])
+            @test p - p == P([0*a])
+
+            # poly mult
+            @test p * q == P(conv([a,b,c], [a,b]))
+            @test q * p == P(conv([a,b], [a,b, c]))
+
+            # poly powers
+            @test p^2 == p * p
+
+            # evalution
+            @test p(s) == a + b * s + c * s * s
+            @test p(c) == a + b * c + c * c * c
+
+            # ∂, ∫
+            @test derivative(p) == P([b, 2c])
+            @test integrate(p) == P([0*a, a, b/2, c/3])
+
+            # matrix element
+            pq = [p q]
+            @test pq[1] == p
+            @test pq[2] == q
+
+            # implicit promotion
+            @test p + s == P([a+s, b, c])
+            @test_throws Union{ArgumentError, MethodError} p + d == P([a+d, b, c]) # can't fix
+            @test p + P([d]) == P([a+d,b,c])
+
+            ps = [p s]
+            @test ps[1] == p
+            @test ps[2] == s
+        end
+    end
+
+    @testset "T=Matrix (2x2)" begin
+        for P ∈ (Polynomial, ImmutablePolynomial)
+            a,b,c = [1 0; 1 1], [1 0; 2 1], [1 0; 3 1]
+            p = P([a,b,c])
+            q = P([a,b])
+            s = 2
+            d = [4 1; 1 0]
+
+            # scalar product
+            @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]])
+            @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]])
+            @test d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]])
+            @test p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]])
+
+            # poly add
+            @test p + q == P([a+a,b+b,c])
+            @test p - q == P([a-a,b-b,c])
+            @test_throws MethodError p - p == P([0*a])  # no zeros to make zero polynomial
+
+            # poly mult
+            @test p * q == P(conv([a,b,c], [a,b]))
+            @test q * p == P(conv([a,b], [a,b, c]))
+
+            # poly powers
+            @test p^2 == p * p
+
+            # evalution
+            @test p(s) == a + b * s + c * s * s
+            @test p(c) == a + b * c + c * c * c
+
+            # ∂, ∫
+            @test derivative(p) == P([b, 2c])
+            @test integrate(p) == P([0*a, a, b/2, c/3])
+
+            # matrix element
+            @test [p q][1] == p
+            @test [p q][2] == q
+
+            # implicit promotion
+            @test_throws MethodError p + s == P([a+s, b, c]) # OK, no a + s
+            @test p + d == P([a+d, b, c])
+            @test p + P([d]) == P([a+d,b,c])
+
+            @test_throws MethodError [p s][1] == p # no promotion T(s)
+            @test_throws MethodError [p s][2] == s
+        end
+    end
+
+
+    @testset "T=Vector{Int}" begin
+        for P ∈ (Polynomial, ImmutablePolynomial)
+            a,b,c = [1,0,0], [1,1,0], [1,1,1]
+            p = P([a,b,c])
+            q = P([a,b])
+            s = 2
+            d = [1,2,3]
+
+            # scalar product
+            @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]])
+            @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]])
+            @test_throws MethodError d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]]) # Ok, no * for T
+            @test_throws MethodError p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]]) # Ok, no * for T
+
+            # poly add
+            @test p + q == P([a+a,b+b,c])
+            @test p - q == P([a-a,b-b,c])
+            @test_throws MethodError p - p == P([0*a])  # no zero(T) to make zero polynomial
+
+            # poly mult
+            @test_throws MethodError p * q == P(conv([a,b,c], [a,b])) # Ok, no * for T
+            @test_throws MethodError q * p == P(conv([a,b], [a,b, c])) # Ok, no * for T
+
+            # poly powers
+            @test_throws MethodError p^2 == p * p # Ok, no * for T
+
+            # evalution
+            @test p(s) == a + b * s + c * s * s
+            @test_throws MethodError p(c) == a + b * c + c * c * c # OK, no b * c
+
+            # ∂, ∫
+            @test derivative(p) == P([b, 2c])
+            @test integrate(p) == P([0*a, a, b/2, c/3])
+
+
+            # matrix element
+            @test [p q][1] == p
+            @test [p q][2] == q
+
+            # implicit promotion
+            @test_throws MethodError p + s == P([a+s, b, c])  # OK, no a + s
+            @test  p + d == P([a+d, b, c])
+            @test p + P([d]) == P([a+d,b,c])
+
+            @test_throws MethodError [p s][1] == p # no promotion T(s)
+            @test_throws MethodError [p s][2] == s
+        end
+    end
+
+
+    if VERSION >= v"1.5"
+        eval(quote
+             using StaticArrays
+             end)
+        @testset "T=SA" begin
+            for P ∈ (Polynomial, ImmutablePolynomial )
+                a,b,c = SA[1 0; 1 1], SA[1 0; 2 1], SA[1 0; 3 1]
+                p = P([a,b,c])
+                q = P([a,b])
+                s = 2
+                d = SA[4 1; 1 0]
+
+                # scalar product
+                @test s*p == P([s*cᵢ for cᵢ ∈ [a,b,c]])
+                @test p*s == P([cᵢ*s for cᵢ ∈ [a,b,c]])
+                @test d*p == P([d*cᵢ for cᵢ ∈ [a,b,c]])
+                @test p*d == P([cᵢ*d for cᵢ ∈ [a,b,c]])
+
+                # poly add
+                @test p + q == P([a+a,b+b,c])
+                @test p - p == P([0*a])
+
+                # poly mult
+                @test p * q == P(conv([a,b,c], [a,b]))
+                @test q * p == P(conv([a,b], [a,b, c]))
+
+                # poly powers
+                @test p^2 == p * p
+
+
+                # evalution
+                @test p(s) == a + b * s + c * s * s
+                @test p(c) == a + b * c + c * c * c
+
+                # ∂, ∫
+                @test derivative(p) == P([b, 2c])
+                @test integrate(p) == P([0*a, a, b/2, c/3])
+
+                # matrix element
+                @test [p q][1] == p
+                @test [p q][2] == q
+
+                # implicit promotion
+                # @test_broken p + s == P([a .+ s, b, c]) # should error, doesn't
+                @test p + d == P([a + d, b, c])
+                @test p + P([d]) == P([a + d,b,c])
+
+                @test_throws MethodError [p s][1] == p # no promotion T(s)
+                @test_throws MethodError [p s][2] == s #
+            end
+        end
+    end
+end
+end
+
 @testset "OffsetVector" begin
     as = ones(3:4)
     bs = parent(as)
-    
-    
+
+
     for P in Ps
         # LaurentPolynomial accepts OffsetArrays; others throw warning
         if P == LaurentPolynomial
@@ -154,11 +362,11 @@ end
             @test P(as) == P(bs)
             @test P{eltype(as)}(as) == P{eltype(as)}(bs)
             # (Or throw an error?)
-            # @test_throws ArgumentError P(as) 
+            # @test_throws ArgumentError P(as)
             # @test P{eltype(as)}(as) == P{eltype(as)}(bs)
         end
     end
-        
+
     a = [1,1]
     b = OffsetVector(a, axes(a))
     c = ZVector(a)
@@ -167,7 +375,7 @@ end
         if P == LaurentPolynomial && continue
             @test P(a) == P(b) == P(c) == P(d)
         end
-        
+
     end
 end
 
@@ -266,7 +474,7 @@ end
         @test pX != pS1
         @test pS1 == pS2
         @test pS1 == pS3
-        
+
         @test indeterminate(pS1 + pS1) == indeterminate(pS1)
         @test indeterminate(pS1 - pS1) == indeterminate(pS1)
         @test indeterminate(pS1 * pS1) == indeterminate(pS1)
@@ -446,8 +654,8 @@ end
         @test Polynomials.evalpoly.(1/2, ps) ≈ [p(1/2)  for  p  in ps]
 
     end
-    
-    
+
+
     # constant polynomials and type
     Ts = (Int, Float32, Float64, Complex{Int}, Complex{Float64})
     for P in (Polynomial, ImmutablePolynomial, SparsePolynomial)
@@ -761,21 +969,21 @@ end
     end
     meths = (Base.vect, Base.vcat, Base.hcat)
     for P in (Polynomial, ImmutablePolynomial, SparsePolynomial, LaurentPolynomial)
-        
+
         p,q = P([1,2], :x), P([1,2], :y)
         P′′ = P == LaurentPolynomial ? P : P′ # different promotion rule
-        
+
         # * should promote to Polynomial type if mixed (save Laurent Polynomial)
         @testset "promote mixed polys" begin
             for m ∈ meths
-                @test _test(m(p,p), P, :x)            
+                @test _test(m(p,p), P, :x)
                 @test _test(m(p,r), P′′, :x)
             end
-            
+
             @test _test(Base.hvcat((2,1), p, r,[p r]), P′′, :x)
-            
+
         end
-        
+
         # * numeric constants should promote to a polynomial, when mixed
         @testset "promote numbers to poly" begin
             for m ∈ meths
@@ -784,20 +992,20 @@ end
                 @test _test(m(1,1,p), P, :x)
                 @test _test(m(p,1,1), P, :x)
             end
-            
-            @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x)        
+
+            @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x)
         end
-        
+
         # * non-constant polynomials must share the same indeterminate
         @testset "non constant polys share same X" begin
             for m ∈ meths
                 @test_throws ArgumentError m(p,q)
                 @test_throws ArgumentError m(p,s)
             end
-            
+
             @test_throws ArgumentError Base.hvcat((2,1), p, q,[p q])
         end
-        
+
 
         # * constant polynomials are treated as `P{T,X}`, not elements of `T`
          @testset "constant polys" begin
@@ -807,10 +1015,10 @@ end
                 @test _test(m(1,1,one(p)), P, :x)
                 @test _test(m(one(p),1,1), P, :x)
             end
-            
-            @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x)        
+
+            @test _test(Base.hvcat((3,1), 1, p, r,[1 p r]), P′′, :x)
         end
-        
+
         # * Promotion can be forced to mix constant-polynomials
         @testset "Use typed constructor to mix constant polynomals" begin
             𝑷,𝑸 = P{Int,:x}, P{Int,:y} # not typeof(p),... as Immutable carries N
@@ -827,10 +1035,10 @@ end
             @test eltype(𝑷[1 one(p); one(p) one(q)]) == 𝑷
             @test eltype(𝑸[1 one(p); one(p) one(q)]) == 𝑸
         end
-        
+
         @testset "hvcat" begin
             p,q = P([1,2],:x), P([1,2],:y)
-            
+
             q1 = [q 1]
             q11 = [one(q) 1]
 
@@ -843,10 +1051,10 @@ end
             @test eltype(𝑷[1 p; q11]) == 𝑷
             @test eltype(Base.typed_hvcat(𝑷, (2, 1), 1, p, q11)) == 𝑷
         end
-        
+
     end
 
-    
+
 end
 
 @testset "Linear Algebra" begin
@@ -1143,7 +1351,7 @@ end
     p = Polynomial([1.234567890, 2.34567890])
     io=IOBuffer(); printpoly(io, p, compact=true); @test String(take!(io)) == "1.23457 + 2.34568*x"
     io=IOBuffer(); printpoly(io, p, compact=true, mulsymbol=""); @test String(take!(io)) == "1.23457 + 2.34568x"
-    
+
     ## issue 278 with complex
     @test printpoly_to_string(Polynomial([1 + im, 1, 2, im, 2im, 1+im, 1-im])) == "1 + im + x + 2*x^2 + im*x^3 + 2im*x^4 + (1 + im)x^5 + (1 - im)x^6"