Eltype fix; use T[] for 0-polynomial; speed up LaurentPolynomial operations (#487)

* T[] for zero polynomial
* trim out use of one(T)
* speed up Laurent;

Author: jverzani

docs/src/reference.md

@@ -79,7 +79,7 @@ vander
 
 ## Plotting
 
-Polynomials can be plotted directly using [Plots.jl](https://github.com/juliaplots/plots.jl).
+Polynomials can be plotted directly using [Plots.jl](https://github.com/juliaplots/plots.jl) or [Makie.jl](https://github.com/MakieOrg/Makie.jl).
 
 ```julia
 plot(::AbstractPolynomial; kwds...)

src/abstract.jl

@@ -30,8 +30,8 @@ Some `T`s will not be successful
 * 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})`
+* add/sub: `p::Polynomial{T} + q::Polynomial{T}` should be defined
+* sub: `p -p` sometimes needs `zero{T}` defined
 
 * 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
@@ -41,6 +41,7 @@ Some `T`s will not be successful
 
 * evaluation: `p::Polynomial{T}(s::Number)`
 * evaluation `p::Polynomial{T}(c::T)`   needs `T*T` defined
+* evaluation of a `0` polynomial requires `zero(T)` to be defined.
 
 * derivatives: `derivative(p::Polynomial{T})`
 * integrals: `integrate(p::Polynomial{T})`

src/common.jl

@@ -385,7 +385,7 @@ end
 function chop!(ps::Vector{T};
                rtol::Real = Base.rtoldefault(real(T)),
                atol::Real = 0,) where {T}
-
+    isempty(ps) && return ps
     tol = norm(ps) * rtol + atol
     for i = lastindex(ps):-1:1
         val = ps[i]
@@ -705,6 +705,7 @@ degreerange(p::AbstractPolynomial) = firstindex(p):lastindex(p)
 # getindex
 function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T}
     m,M = firstindex(p), lastindex(p)
+    m > M && return zero(T)
     idx < m && throw(BoundsError(p, idx))
     idx > M && return zero(T)
     p.coeffs[idx - m + 1]
@@ -831,7 +832,7 @@ Returns a representation of 0 as the given polynomial.
 """
 function Base.zero(::Type{P}) where {P<:AbstractPolynomial}
     T,X = eltype(P), indeterminate(P)
-    ⟒(P){T,X}(zeros(T,1))
+    ⟒(P){T,X}(T[])
 end
 Base.zero(::Type{P}, var::SymbolLike) where {P <: AbstractPolynomial} = zero(⟒(P){eltype(P),Symbol(var)}) #default 0⋅b₀
 Base.zero(p::P, var=indeterminate(p)) where {P <: AbstractPolynomial} = zero(P, var)
@@ -1031,7 +1032,7 @@ end
 
 ## -- multiplication
 
-
+# Scalar multiplication; no assumption of commutivity
 function scalar_mult(p::P, c::S) where {S, T, X, P<:AbstractPolynomial{T,X}}
     result = coeffs(p) .* (c,)
     ⟒(P){eltype(result), X}(result)
@@ -1044,11 +1045,11 @@ end
 
 scalar_mult(p1::AbstractPolynomial, p2::AbstractPolynomial) = error("scalar_mult(::$(typeof(p1)), ::$(typeof(p2))) is not defined.") # avoid ambiguity, issue #435
 
-function Base.:/(p::P, c::S) where {P <: AbstractPolynomial,S}
-    _convert(p, coeffs(p) ./ c)
-end
+# scalar div
+Base.:/(p::P, c::S) where {P <: AbstractPolynomial,S} = scalar_div(p, c)
+scalar_div(p::AbstractPolynomial, c) = scalar_mult(p, inv(c))
 
-## polynomial p*q
+## Polynomial p*q
 ## Polynomial multiplication formula depend on the particular basis used. The subtype must implement
 function Base.:*(p1::P, p2::Q) where {T,X,P <: AbstractPolynomial{T,X},S,Y,Q <: AbstractPolynomial{S,Y}}
     isconstant(p1) && return constantterm(p1) * p2
@@ -1060,6 +1061,7 @@ end
 
 Base.:^(p::AbstractPolynomial, n::Integer) = Base.power_by_squaring(p, n)
 
+
 function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial}
     n, d = promote(num, den)
     return divrem(n, d)

src/contrib.jl

@@ -5,6 +5,8 @@ using Base.Cartesian
 
 # direct version (do not check if threshold is satisfied)
 function fastconv(E::Array{T,N}, k::Array{T,N}) where {T,N}
+    isempty(E) && return E
+    isempty(k) && return k
     retsize = ntuple(n -> size(E, n) + size(k, n) - 1, Val{N}())
     ret = zeros(T, retsize)
     convn!(ret, E, k)
@@ -34,7 +36,6 @@ end
 module EvalPoly
 using LinearAlgebra
 function evalpoly(x::S, p::Tuple) where {S}
-    p == () && return zero(S)
     if @generated
         N = length(p.parameters)
         ex = :(p[end]*_one(S))
@@ -51,14 +52,13 @@ evalpoly(x, p::AbstractVector) = _evalpoly(x, p)
 
 # https://discourse.julialang.org/t/i-have-a-much-faster-version-of-evalpoly-why-is-it-faster/79899; improvement *and* closes #313
 function _evalpoly(x::S, p) where {S}
-    i = lastindex(p)
 
-    @inbounds out = p[i]*_one(x)
+    i = lastindex(p)
+    @inbounds out = p[i] * _one(x)
     i -= 1
-
     while i >= firstindex(p)
 	@inbounds out = _muladd(out, x, p[i])
-		i -= 1
+	i -= 1
     end
 
     return out
@@ -68,8 +68,8 @@ end
 function evalpoly(z::Complex, p::Tuple)
     if @generated
         N = length(p.parameters)
-        a = :(p[end]*_one(z))
-        b = :(p[end-1]*_one(z))
+        a = :(p[end] .+ _zero(z))  # avoid one(x)
+        b = :(p[end-1] .+ _zero(z))
         as = []
         for i in N-2:-1:1
             ai = Symbol("a", i)
@@ -124,6 +124,14 @@ _muladd(a, b::Matrix, c) = (a*I)*b + c*I
 _one(P::Type{<:Matrix}) = one(eltype(P))*I
 _one(x::Matrix) = one(eltype(x))*I
 _one(x) = one(x)
+
+_zero(P::Type{<:Matrix}) = zero(eltype(P))*I
+function _zero(x::Matrix)
+    m = LinearAlgebra.checksquare(x)
+    zero(eltype(x)) * I(m)
+end
+_zero(x) = zero(x)
+
 end
 
 ## get type of parametric composite type without type parameters

src/polynomials/ImmutablePolynomial.jl

@@ -64,7 +64,7 @@ end
 
 ## Various interfaces
 ## Abstract Vector coefficients
-function ImmutablePolynomial{T,X, N}(coeffs::AbstractVector{T})  where {T,X, N}
+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)
 
@@ -188,44 +188,36 @@ function Base.:+(p1::P, p2::Q) where {T,X,N,P<:ImmutablePolynomial{T,X,N},
 end
 
 ## multiplication
-
 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
-    return ImmutablePolynomial(cs, X)
- end
+    iszero(N) && return zero(ImmutablePolynomial{promote_type(T,S),X})
+    iszero(c) && ImmutablePolynomial([p[0] .* c], X)
+    return _polynomial(p, p.coeffs .* (c,))
+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)
+    iszero(N) && return zero(ImmutablePolynomial{promote_type(T,S),X})
+    iszero(c) && ImmutablePolynomial([c .* p[0]], X)
+    return _polynomial(p, (c,) .* p.coeffs)
 end
 
-function Base.:/(p::ImmutablePolynomial{T,X,N}, c::S) where {T,X,N,S<:Number}
-    R = eltype(one(T)/one(S))
+function _polynomial(p::ImmutablePolynomial{T,X,N}, cs) where {T, X, N}
+    R = eltype(cs)
     P = ImmutablePolynomial{R,X}
-    (N == 0  || isinf(c)) && return zero(P)
-    cs = p.coeffs ./ c
-    iszero(cs[end]) ? P(cs) : P{N}(cs) # more performant to specify when N is known
+    iszero(cs[end]) ? P(cs) : P{N}(cs)
 end
 
-
 function Base.:*(p1::ImmutablePolynomial{T,X,N}, p2::ImmutablePolynomial{S,X,M}) where {T,S,X,N,M}
-    R = promote_type(T,S)
-    P = ImmutablePolynomial{R,X}
-
-    (iszero(N) || iszero(M)) && return zero(P)
+    (iszero(N) || iszero(M)) && return zero(ImmutablePolynomial{promote_type(T,S),X})
 
     cs = ⊗(ImmutablePolynomial, p1.coeffs, p2.coeffs) #(p1.coeffs) ⊗ (p2.coeffs)
+    R = eltype(cs)
+    P = ImmutablePolynomial{R,X}
     iszero(cs[end]) ? P(cs) : P{N+M-1}(cs)  # more performant to specify when N is known
 
 end
 
 Base.to_power_type(p::ImmutablePolynomial{T,X,N}) where {T,X,N} = p
 
-
 ## more performant versions borrowed from StaticArrays
 ## https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/linalg.jl
 LinearAlgebra.norm(q::ImmutablePolynomial{T,X,0}) where {T,X} = zero(real(float(T)))

src/polynomials/LaurentPolynomial.jl

@@ -1,5 +1,6 @@
 export LaurentPolynomial
 
+
 """
     LaurentPolynomial{T,X}(coeffs::AbstractVector, [m::Integer = 0], [var = :x])
 
@@ -94,6 +95,11 @@ struct LaurentPolynomial{T, X} <: LaurentBasisPolynomial{T, X}
             new{T,X}(c, Ref(m′),  Ref(n))
         end
     end
+    # non copying version assumes trimmed coeffs
+    function LaurentPolynomial{T,X}(::Val{false}, coeffs::AbstractVector{T},
+                                    m::Integer=0) where {T, X}
+        new{T,X}(coeffs, Ref(m), Ref(m + length(coeffs) - 1))
+    end
 end
 
 @register LaurentPolynomial
@@ -259,7 +265,7 @@ end
 function showterm(io::IO, ::Type{<:LaurentPolynomial}, pj::T, var, j, first::Bool, mimetype) where {T}
     if iszero(pj) return false end
     pj = printsign(io, pj, first, mimetype)
-    if !(pj == one(T) && !(showone(T) || j == 0))
+    if !(hasone(T) && pj == one(T) && !(showone(T) || j == 0))
         printcoefficient(io, pj, j, mimetype)
     end
     printproductsign(io, pj, j, mimetype)
@@ -413,57 +419,114 @@ function Base.:+(p::LaurentPolynomial{T,X}, c::S) where {T, X, S <: Number}
 end
 
 ##
-## Poly + and  *
-##
-function Base.:+(p1::P, p2::P) where {T,X,P<:LaurentPolynomial{T,X}}
+## Poly +, - and  *
+## uses some ideas from https://github.com/jmichel7/LaurentPolynomials.jl/blob/main/src/LaurentPolynomials.jl for speedups
+Base.:+(p1::LaurentPolynomial, p2::LaurentPolynomial) = add_sub(+, p1, p2)
+Base.:-(p1::LaurentPolynomial, p2::LaurentPolynomial) = add_sub(-, p1, p2)
 
-    isconstant(p1) && return constantterm(p1) + p2
-    isconstant(p2) && return p1 + constantterm(p2)
+function add_sub(op, p1::P, p2::Q) where {T, X, P <: LaurentPolynomial{T,X},
+                                          S, Y, Q <: LaurentPolynomial{S,Y}}
 
+    isconstant(p1) && return op(constantterm(p1), p2)
+    isconstant(p2) && return op(p1, constantterm(p2))
+    assert_same_variable(X, Y)
 
     m1,n1 = (extrema ∘ degreerange)(p1)
     m2,n2 = (extrema ∘ degreerange)(p2)
-    m,n = min(m1,m2), max(n1, n2)
+    m, n = min(m1,m2), max(n1, n2)
 
-    as = zeros(T, length(m:n))
-    for i in m:n
-        as[1 + i-m] = p1[i] + p2[i]
-    end
+    R = typeof(p1.coeffs[1] + p2.coeffs[1]) # non-empty
+    as = zeros(R, length(m:n))
+
+    d = m1 - m2
+    d1, d2 = m1 > m2 ? (d,0) : (0, -d)
 
-    q = P(as, m)
-    chop!(q)
+    for (i, pᵢ) ∈ pairs(p1.coeffs)
+        @inbounds as[d1 + i] = pᵢ
+    end
+    for (i, pᵢ) ∈ pairs(p2.coeffs)
+        @inbounds as[d2 + i] = op(as[d2+i], pᵢ)
+    end
 
+    m = _laurent_chop!(as, m)
+    isempty(as) && return zero(LaurentPolynomial{R,X})
+    q = LaurentPolynomial{R,X}(Val(false), as, m)
     return q
 
 end
 
-function Base.:*(p1::P, p2::P) where {T,X,P<:LaurentPolynomial{T,X}}
+function Base.:*(p1::P, p2::Q) where {T,X,P<:LaurentPolynomial{T,X},
+                                      S,Y,Q<:LaurentPolynomial{S,Y}}
+
+    isconstant(p1) && return constantterm(p1) * p2
+    isconstant(p2) && return p1 * constantterm(p2)
+    assert_same_variable(X, Y)
 
     m1,n1 = (extrema ∘ degreerange)(p1)
     m2,n2 = (extrema ∘ degreerange)(p2)
     m,n = m1 + m2, n1+n2
 
-    as = zeros(T, length(m:n))
-    for i in eachindex(p1)
-        p1ᵢ = p1[i]
-        for j in eachindex(p2)
-            as[1 + i+j - m] = muladd(p1ᵢ, p2[j], as[1 + i + j - m])
+    R = promote_type(T,S)
+    as = zeros(R, length(m:n))
+
+    for (i, p₁ᵢ) ∈ pairs(p1.coeffs)
+        for (j, p₂ⱼ) ∈ pairs(p2.coeffs)
+            @inbounds as[i+j-1] += p₁ᵢ * p₂ⱼ
         end
     end
 
-    p = P(as, m)
-    chop!(p)
+    m = _laurent_chop!(as, m)
+
+    isempty(as) && return zero(LaurentPolynomial{R,X})
+    p = LaurentPolynomial{R,X}(Val(false), as, m)
 
     return p
 end
 
+function _laurent_chop!(as, m)
+    while !isempty(as)
+        if iszero(first(as))
+            m += 1
+            popfirst!(as)
+        else
+            break
+        end
+    end
+    while !isempty(as)
+        if iszero(last(as))
+            pop!(as)
+        else
+            break
+        end
+    end
+    m
+end
+
 function scalar_mult(p::LaurentPolynomial{T,X}, c::Number) where {T,X}
     LaurentPolynomial(p.coeffs .* c, p.m[], Var(X))
 end
 function scalar_mult(c::Number, p::LaurentPolynomial{T,X}) where {T,X}
     LaurentPolynomial(c .* p.coeffs, p.m[], Var(X))
 end
 
+
+function integrate(p::P) where {T, X, P <: LaurentBasisPolynomial{T, X}}
+
+    R = typeof(constantterm(p) / 1)
+    Q = ⟒(P){R,X}
+
+    hasnan(p) && return Q([NaN])
+    iszero(p) && return zero(Q)
+
+    ∫p = zero(Q)
+    for (k, pₖ) ∈ pairs(p)
+        iszero(pₖ) && continue
+        k == -1 && throw(ArgumentError("Can't integrate Laurent polynomial with  `x⁻¹` term"))
+        ∫p[k+1] = pₖ/(k+1)
+    end
+    ∫p
+end
+
 ##
 ## roots
 ##

src/polynomials/Poly.jl

@@ -116,7 +116,7 @@ end
 
 function Polynomials.integrate(p::P, k::S) where {T, X, P <: Poly{T, X}, S<:Number}
 
-    R = eltype((one(T)+one(S))/1)
+    R = eltype(one(T)/1 + one(S))
     Q = Poly{R,X}
     if hasnan(p) || isnan(k)
         return P([NaN]) # keep for Poly, not Q