Negative exponents in SparsePolynomial (#411)

* negative exponents in SparsePolynomials

* fix conversions to polynomial types

* fix conversions

* version bump to v3.1.2

* test conversions

* test minimumexponent for SparsePolynomial

* test show for negative exponents

* firstindex instead of minimumexponent in LP

* missing tests for PolynomialValues/Keys

* fix testset

* non-allocating == check

Author: jishnub

Project.toml

@@ -2,7 +2,7 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "3.1.1"
+version = "3.1.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

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,9 +612,15 @@ 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)
 
 # getindex
@@ -671,7 +676,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 +685,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,8 +1024,15 @@ Base.rem(n::AbstractPolynomial, d::AbstractPolynomial) = divrem(n, d)[2]
 #=
 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))
+function Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial)
+    check_same_variable(p1,p2) || return false
+    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
 Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && constantterm(p) == n
 Base.:(==)(n::Number, p::AbstractPolynomial) = p == n
 

src/polynomials/LaurentPolynomial.jl

@@ -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,9 +209,10 @@ 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)
+minimumexponent(::Type{<:LaurentPolynomial}) = typemin(Int)
+minimumexponent(p::LaurentPolynomial) = p.m[]
+Base.firstindex(p::LaurentPolynomial) = minimumexponent(p)
+degree(p::LaurentPolynomial) = p.n[]
 degreerange(p::LaurentPolynomial) = firstindex(p):lastindex(p)
 
 _convert(p::P, as) where {T,X,P <: LaurentPolynomial{T,X}} = ⟒(P)(as, firstindex(p), Var(X))

src/polynomials/SparsePolynomial.jl

@@ -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
 
@@ -141,7 +135,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 +195,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 +224,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)

src/polynomials/standard-basis.jl

@@ -49,9 +49,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
 

test/Poly.jl

@@ -1,4 +1,13 @@
 ## Tests for compatiability with the old form
+# Create a separate module so that include(file) doesn't import PolyCompat to Main
+module PolyTests
+
+using LinearAlgebra
+using Test
+using SparseArrays
+using ..Polynomials
+using ..SpecialFunctions
+
 # test for not loaded
 @test_throws UndefVarError Poly([1,2,3])
 @test_throws UndefVarError poly([1,2,3])
@@ -29,7 +38,7 @@ p = Polynomial([1,2,3])
 @test_throws ErrorException polyder(p)
 
 
-        
+
 
 
 ## ---------
@@ -280,15 +289,15 @@ repr(p)
 
 ## changes to show
 p = Poly([1,2,3,1])  # leading coefficient of 1
-@test repr(p) == "Poly(1 + 2*x + 3*x^2 + x^3)"
+@test repr(p) == "$Poly(1 + 2*x + 3*x^2 + x^3)"
 p = Poly([1.0, 2.0, 3.0, 1.0])
-@test repr(p) == "Poly(1.0 + 2.0*x + 3.0*x^2 + 1.0*x^3)"
+@test repr(p) == "$Poly(1.0 + 2.0*x + 3.0*x^2 + 1.0*x^3)"
 p = Poly([1, im])
-@test repr(p) == "Poly(1 + im*x)"
+@test repr(p) == "$Poly(1 + im*x)"
 p = Poly([1+im, 1-im, -1+im, -1 - im])# minus signs
-@test repr(p) == "Poly(1 + im + (1 - im)x - (1 - im)x^2 - (1 + im)x^3)"
+@test repr(p) == "$Poly(1 + im + (1 - im)x - (1 - im)x^2 - (1 + im)x^3)"
 p = Poly([1.0, 0 + NaN*im, NaN, Inf, 0 - Inf*im]) # handle NaN or Inf appropriately
-@test repr(p) == "Poly(1.0 + NaN*im*x + NaN*x^2 + Inf*x^3 - Inf*im*x^4)"
+@test repr(p) == "$Poly(1.0 + NaN*im*x + NaN*x^2 + Inf*x^3 - Inf*im*x^4)"
 
 p = Poly([1,2,3])
 @test repr("text/latex", p) == "\$1 + 2\\cdot x + 3\\cdot x^{2}\$"
@@ -497,3 +506,5 @@ end
     PQexpint = Pade(d, 30, 30)
     @test Float64(PQexpint(1.0)) ≈ exp(1) * (-γ - sum([(-1)^k / k / gamma(k + 1) for k = 1:20]))
 end
+
+end # module