fix == and isapprox

Author: jverzani

src/common.jl

@@ -407,7 +407,6 @@ Calculates the p-norm of the polynomial's coefficients
 """
 function LinearAlgebra.norm(q::AbstractPolynomial, p::Real = 2)
     vs = values(q)
-    isempty(vs) && return zero(eltype(q))
     return norm(vs, p) # if vs=() must be handled in special type
 end
 
@@ -572,7 +571,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(coeffs(p)) ? -1 : length(coeffs(p)) - 1 + min(0, minimumexponent(p))
 
 
 """
@@ -613,10 +612,17 @@ mapdomain(::P, x::AbstractArray) where {P <: AbstractPolynomial} = mapdomain(P,
 
 #=
 indexing =#
+# minimumexponent(p) returns min(0, minimum(keys(p)))
+# For most polynomial types, this is statically known to be zero
+# For polynomials that support negative indices, minimumexponent(typeof(p))
+# should return typemin(Int)
+minimumexponent(p::AbstractPolynomial) = minimumexponent(typeof(p))
+minimumexponent(::Type{<:AbstractPolynomial}) = 0
 Base.firstindex(p::AbstractPolynomial) = 0
-Base.lastindex(p::AbstractPolynomial) = length(p) - 1
-Base.eachindex(p::AbstractPolynomial) = 0:length(p) - 1
+Base.lastindex(p::AbstractPolynomial) = length(p) - 1 + firstindex(p)
+Base.eachindex(p::AbstractPolynomial) = firstindex(p):lastindex(p)
 Base.broadcastable(p::AbstractPolynomial) = Ref(p)
+degreerange(p::AbstractPolynomial) = firstindex(p):lastindex(p)
 
 # getindex
 function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T}
@@ -671,7 +677,7 @@ Base.iterate(p::AbstractPolynomial, state = firstindex(p)) = _iterate(p, state)
 struct PolynomialKeys{P} <: AbstractSet{Int}
     p::P
 end
-struct PolynomialValues{P, T} <: AbstractSet{T}
+struct PolynomialValues{P, T}
     p::P
 
     PolynomialValues{P}(p::P) where {P} = new{P, eltype(p)}(p)
@@ -680,6 +686,7 @@ end
 Base.keys(p::AbstractPolynomial) =  PolynomialKeys(p)
 Base.values(p::AbstractPolynomial) =  PolynomialValues(p)
 Base.length(p::PolynomialValues) = length(p.p.coeffs)
+Base.eltype(p::PolynomialValues{<:Any,T}) where {T} = T
 Base.length(p::PolynomialKeys) = length(p.p.coeffs)
 Base.size(p::Union{PolynomialValues, PolynomialKeys}) = (length(p),)
 function Base.iterate(v::PolynomialKeys, state = firstindex(v.p))
@@ -1018,23 +1025,33 @@ Base.rem(n::AbstractPolynomial, d::AbstractPolynomial) = divrem(n, d)[2]
 #=
 Comparisons =#
 Base.isequal(p1::P, p2::P) where {P <: AbstractPolynomial} = hash(p1) == hash(p2)
+function Base.:(==)(p1::P, p2::P) where {P <: AbstractPolynomial}
+    iszero(p1) && iszero(p2) && return true
+    eachindex(p1) == eachindex(p2) || return false
+    # coeffs(p1) == coeffs(p2), but non-allocating
+    p1val = (p1[i] for i in eachindex(p1))
+    p2val = (p2[i] for i in eachindex(p2))
+    all(((a,b),) -> a == b, zip(p1val, p2val))
+end
 function Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial)
     if isconstant(p1)
         isconstant(p2) && return constantterm(p1) == constantterm(p2)
         return false
+    elseif isconstant(p2)
+        return false # p1 is not constant
     end
-    indeterminate(p1) == indeterminate(p2) || return false
+    check_same_variable(p1, p2) || return false
     ==(promote(p1,p2)...)
 end
-Base.:(==)(p1::P, p2::P) where {P <: AbstractPolynomial} =
-    check_same_variable(p1,p2) && (coeffs(p1) == coeffs(p2))
 Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && constantterm(p) == n
 Base.:(==)(n::Number, p::AbstractPolynomial) = p == n
 
 function Base.isapprox(p1::AbstractPolynomial, p2::AbstractPolynomial; kwargs...)
     if isconstant(p1)
         isconstant(p2) && return constantterm(p1) == constantterm(p2)
         return false
+    elseif isconstant(p2)
+        return false
     end
     isapprox(promote(p1, p2)...; kwargs...)
 end

src/polynomials/LaurentPolynomial.jl

@@ -51,7 +51,7 @@ julia> integrate(q)
 LaurentPolynomial(1.0*x + 0.5*x² + 0.3333333333333333*x³)
 
 julia> integrate(p)  # x⁻¹  term is an issue
-ERROR: ArgumentError: Can't integrate Laurent  polynomial with  `x⁻¹` term
+ERROR: ArgumentError: Can't integrate Laurent polynomial with  `x⁻¹` term
 
 julia> integrate(P([1,1,1], -5))
 LaurentPolynomial(-0.25*x⁻⁴ - 0.3333333333333333*x⁻³ - 0.5*x⁻²)
@@ -69,7 +69,7 @@ julia> x^degree(p) * p(x⁻¹) # reverses  coefficients
 LaurentPolynomial(3.0 + 2.0*x + 1.0*x²)
 ```
 """
-struct LaurentPolynomial{T, X} <: StandardBasisPolynomial{T, X}
+struct LaurentPolynomial{T, X} <: LaurentBasisPolynomial{T, X}
     coeffs::Vector{T}
     m::Base.RefValue{Int}
     n::Base.RefValue{Int}
@@ -126,12 +126,6 @@ Base.promote_rule(::Type{P},::Type{Q}) where {T, X, P <: LaurentPolynomial{T,X},
 Base.promote_rule(::Type{Q},::Type{P}) where {T, X, P <: LaurentPolynomial{T,X}, S, Q <: StandardBasisPolynomial{S,X}} =
     LaurentPolynomial{promote_type(T, S),X}
 
-function Base.convert(P::Type{<:Polynomial}, q::LaurentPolynomial)
-    m,n = (extrema∘degreerange)(q)
-    m < 0 && throw(ArgumentError("Can't convert a Laurent polynomial with m < 0"))
-    P([q[i] for i  in 0:n], indeterminate(q))
-end
-
 # need to add p.m[], so abstract.jl method isn't sufficent
 # XXX unlike abstract.jl, this uses Y variable in conversion; no error
 # Used in DSP.jl
@@ -154,7 +148,7 @@ function Base.convert(::Type{P}, q::StandardBasisPolynomial{S}) where {P <:Laure
 
      T = _eltype(P, q)
      X = indeterminate(P, q)
-     ⟒(P){T,X}([q[i] for i in 0:degree(q)], 0)
+     ⟒(P){T,X}([q[i] for i in eachindex(q)], firstindex(q))
  end
 
 function Base.convert(::Type{P}, q::AbstractPolynomial) where {P <:LaurentPolynomial}
@@ -215,10 +209,11 @@ function Base.setindex!(p::LaurentPolynomial{T}, value::Number, idx::Int) where
 
 end
 
-Base.firstindex(p::LaurentPolynomial) = p.m[]
-Base.lastindex(p::LaurentPolynomial) = p.n[]
-Base.eachindex(p::LaurentPolynomial) = degreerange(p)
-degreerange(p::LaurentPolynomial) = firstindex(p):lastindex(p)
+minimumexponent(::Type{<:LaurentPolynomial}) = typemin(Int)
+minimumexponent(p::LaurentPolynomial) = p.m[]
+Base.firstindex(p::LaurentPolynomial) = minimumexponent(p)
+degree(p::LaurentPolynomial) = p.n[]
+
 
 _convert(p::P, as) where {T,X,P <: LaurentPolynomial{T,X}} = ⟒(P)(as, firstindex(p), Var(X))
 
@@ -517,12 +512,13 @@ function derivative(p::P, order::Integer = 1) where {T, X, P<:LaurentPolynomial{
     n = n - order
     as =  zeros(T, length(m:n))
 
-    for k in eachindex(p)
+    for (k, pₖ) in pairs(p)
+        iszero(pₖ) && continue
         idx = 1 + k - order - m
         if 0 ≤ k ≤ order - 1
             as[idx] = zero(T)
         else
-            as[idx] = reduce(*, (k - order + 1):k, init = p[k])
+            as[idx] = reduce(*, (k - order + 1):k, init = pₖ)
         end
     end
 
@@ -531,38 +527,6 @@ function derivative(p::P, order::Integer = 1) where {T, X, P<:LaurentPolynomial{
 end
 
 
-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
-
-
-    m,n = (extrema ∘ degreerange)(p)
-    if  n < 0
-        n = 0
-    else
-        n += 1
-    end
-    if m < 0
-        m += 1
-    else
-        m = 0
-    end
-    as = zeros(R,  length(m:n))
-
-    for k in eachindex(p)
-        as[1 + k+1-m]  =  p[k]/(k+1)
-    end
-
-    return Q(as, m)
-
-end
-
 function Base.gcd(p::LaurentPolynomial{T,X}, q::LaurentPolynomial{T,Y}, args...; kwargs...) where {T,X,Y}
     mp, Mp = (extrema ∘ degreerange)(p)
     mq, Mq = (extrema ∘ degreerange)(q)

src/polynomials/SparsePolynomial.jl

@@ -38,7 +38,7 @@ julia> p(1)
 ```
 
 """
-struct SparsePolynomial{T, X} <: StandardBasisPolynomial{T, X}
+struct SparsePolynomial{T, X} <: LaurentBasisPolynomial{T, X}
     coeffs::Dict{Int, T}
     function SparsePolynomial{T, X}(coeffs::AbstractDict{Int, S}) where {T, X, S}
         c = Dict{Int, T}(coeffs)
@@ -73,18 +73,9 @@ function SparsePolynomial{T,X}(coeffs::AbstractVector{S}) where {T, X, S}
     return SparsePolynomial{T,X}(p)
 end
 
-
-
-
-# conversion
-function Base.convert(P::Type{<:Polynomial}, q::SparsePolynomial)
-    ⟒(P)(coeffs(q), indeterminate(q))
-end
-
-function Base.convert(P::Type{<:SparsePolynomial}, q::StandardBasisPolynomial{T}) where {T}
-    R = promote_type(eltype(P), T)
-    ⟒(P){R,indeterminate(P,q)}(coeffs(q))
-end
+minimumexponent(::Type{<:SparsePolynomial}) = typemin(Int)
+minimumexponent(p::SparsePolynomial) = isempty(p.coeffs) ? 0 : min(0, minimum(keys(p.coeffs)))
+Base.firstindex(p::SparsePolynomial) = minimumexponent(p)
 
 ## changes to common
 degree(p::SparsePolynomial) = isempty(p.coeffs) ? -1 : maximum(keys(p.coeffs))
@@ -94,6 +85,8 @@ function isconstant(p::SparsePolynomial)
     return true
 end
 
+Base.convert(::Type{T}, p::StandardBasisPolynomial) where {T<:SparsePolynomial} = T(Dict(pairs(p)))
+
 function basis(P::Type{<:SparsePolynomial}, n::Int)
     T,X = eltype(P), indeterminate(P)
     SparsePolynomial{T,X}(Dict(n=>one(T)))
@@ -104,9 +97,10 @@ end
 function coeffs(p::SparsePolynomial{T})  where {T}
 
     n = degree(p)
-    cs = zeros(T, n+1)
+    cs = zeros(T, length(p))
+    keymin = firstindex(p)
     for (k,v) in p.coeffs
-        cs[k+1]=v
+        cs[k - keymin + 1] = v
     end
     cs
 
@@ -118,7 +112,6 @@ function Base.getindex(p::SparsePolynomial{T}, idx::Int) where {T}
 end
 
 function Base.setindex!(p::SparsePolynomial, value::Number, idx::Int)
-    idx < 0  && return p
     if iszero(value)
         haskey(p.coeffs, idx) && pop!(p.coeffs, idx)
     else
@@ -141,7 +134,10 @@ function Base.iterate(v::PolynomialValues{SparsePolynomial{T,X}}, state...) wher
     return (y[1][2], y[2])
 end
 
-Base.length(S::SparsePolynomial) = isempty(S.coeffs) ? 0 : maximum(keys(S.coeffs)) + 1
+Base.length(S::SparsePolynomial) = isempty(S.coeffs) ? 0 : begin
+    minkey, maxkey = extrema(keys(S.coeffs))
+    maxkey - min(0, minkey) + 1
+end
 
 ##
 ## ----
@@ -198,6 +194,8 @@ function Base.:+(p1::P1, p2::P2) where {T,X, P1<:SparsePolynomial{T,X},
 
 end
 
+Base.:-(a::SparsePolynomial) = typeof(a)(Dict(k=>-v for (k,v) in a.coeffs))
+
 ## Multiplication
 function scalar_mult(p::P, c::S) where {T, X, P <: SparsePolynomial{T,X}, S<:Number}
 
@@ -225,7 +223,6 @@ function scalar_mult(c::S, p::P) where {T, X, P <: SparsePolynomial{T,X}, S<:Num
     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)
@@ -243,32 +240,12 @@ function derivative(p::SparsePolynomial{T,X}, order::Integer = 1) where {T,X}
     hasnan(p) && return P{R,X}(Dict(0 => R(NaN)))
 
     n = degree(p)
-    order > n && return zero(P{R,X})
 
     dpn = zero(P{R,X})
-    @inbounds for k in eachindex(p)
-        dpn[k-order] =  reduce(*, (k - order + 1):k, init = p[k])
+    @inbounds for (k,v) in pairs(p)
+        dpn[k-order] =  reduce(*, (k - order + 1):k, init = v)
     end
 
     return dpn
 
 end
-
-
-function integrate(p::P) where {T, X, P<:SparsePolynomial{T,X}}
-
-    R = eltype(one(T)/1)
-    Q = SparsePolynomial{R,X}
-
-    if hasnan(p)
-        return Q(Dict(0 => NaN))
-    end
-
-    ∫p = zero(Q)
-    for k in eachindex(p)
-        ∫p[k + 1] = p[k] / (k+1)
-    end
-
-    return ∫p
-
-end

src/polynomials/standard-basis.jl

@@ -1,4 +1,5 @@
 abstract type StandardBasisPolynomial{T,X} <: AbstractPolynomial{T,X} end
+abstract type LaurentBasisPolynomial{T,X} <: StandardBasisPolynomial{T,X} end
 
 function showterm(io::IO, ::Type{<:StandardBasisPolynomial}, pj::T, var, j, first::Bool, mimetype) where {T}
 
@@ -49,9 +50,11 @@ function Base.convert(P::Type{<:StandardBasisPolynomial}, q::StandardBasisPolyno
     if isa(q, P)
         return q
     else
+        minimumexponent(P) <= minimumexponent(q) ||
+            throw(ArgumentError("a $P can not have a minimum exponent of $(minimumexponent(q))"))
         T = _eltype(P,q)
         X = indeterminate(P,q)
-        return ⟒(P){T,X}([q[i] for i in 0:degree(q)])
+        return ⟒(P){T,X}([q[i] for i in eachindex(q)])
     end
 end
 
@@ -197,6 +200,22 @@ function integrate(p::P) where {T, X, P <: StandardBasisPolynomial{T, X}}
     return Q(as)
 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
+
 function Base.divrem(num::P, den::Q) where {T, P <: StandardBasisPolynomial{T}, S, Q <: StandardBasisPolynomial{S}}
 
     assert_same_variable(num, den)