Skip to content

Jishnub sparseneg int der #413

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 14 commits into from
May 31, 2022
Merged
1 change: 1 addition & 0 deletions src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down
41 changes: 5 additions & 36 deletions src/polynomials/LaurentPolynomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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))

Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand Down
27 changes: 3 additions & 24 deletions src/polynomials/SparsePolynomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -112,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
Expand Down Expand Up @@ -241,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
17 changes: 17 additions & 0 deletions src/polynomials/standard-basis.jl
Original file line number Diff line number Diff line change
@@ -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}

Expand Down Expand Up @@ -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)
Expand Down
8 changes: 8 additions & 0 deletions test/StandardBasis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down