integrate, derivative fix

Author: jverzani

src/common.jl

@@ -622,6 +622,7 @@ Base.firstindex(p::AbstractPolynomial) = 0
 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}

src/polynomials/LaurentPolynomial.jl

@@ -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}
@@ -213,7 +213,7 @@ 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))
 
@@ -512,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
 
@@ -526,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)
@@ -249,23 +249,3 @@ function derivative(p::SparsePolynomial{T,X}, order::Integer = 1) where {T,X}
     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,v) in pairs(p)
-        k == -1 && throw(ArgumentError("Can't integrate Laurent polynomial with  `x⁻¹` term"))
-        ∫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}
 
@@ -199,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)

test/StandardBasis.jl

@@ -933,6 +933,14 @@ end
         # Issue with overflow and polyder Issue #159
         @test derivative(P(BigInt[0, 1])^100, 100) == P(factorial(big(100)))
     end
+
+    # Sparse and Laurent test issue #409
+    pₛ, pₗ = SparsePolynomial(Dict(-1=>1, 1=>1)), LaurentPolynomial([1,0,1], -1)
+    @test pₛ == pₗ && derivative(pₛ) == derivative(pₗ)
+    @test_throws ArgumentError integrate(pₛ)
+    @test_throws ArgumentError integrate(pₗ)
+    qₛ, qₗ = SparsePolynomial(Dict(-2=>1, 1=>1)), LaurentPolynomial([1,0,0,1], -2)
+    @test qₛ == qₗ && integrate(qₛ) == integrate(qₗ)
 end
 
 @testset "Elementwise Operations" begin