Skip to content

Commit b0adb45

Browse files
blegatblegat
committed
Improvements
1 parent a2faa05 commit b0adb45

File tree

7 files changed

+354
-49
lines changed

7 files changed

+354
-49
lines changed

‎Project.toml‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,10 +5,14 @@ repo = "https://github.com/JuliaAlgebra/MultivariatePolynomials.jl"
55
version = "0.3.2"
66

77
[deps]
8+
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
9+
DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
810
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
911
MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
12+
TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
1013

1114
[compat]
15+
DataStructures = "0.17"
1216
julia = "1"
1317

1418
[extras]

‎src/MultivariatePolynomials.jl‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@ module MultivariatePolynomials
22

33
import LinearAlgebra
44

5+
import DataStructures
6+
57
import MutableArithmetics
68
const MA = MutableArithmetics
79

@@ -13,13 +15,15 @@ Abstract type for a value that can act like a polynomial. For instance, an
1315
`AbstractTerm{T}` is an `AbstractPolynomialLike{T}` since it can act as a
1416
polynomial of only one term.
1517
"""
16-
abstract type AbstractPolynomialLike{T} end
18+
abstract type AbstractPolynomialLike{T} <: MA.AbstractMutable end
19+
1720
"""
1821
AbstractTermLike{T}
1922
2023
Abstract type for a value that can act like a term. For instance, an `AbstractMonomial` is an `AbstractTermLike{Int}` since it can act as a term with coefficient `1`.
2124
"""
2225
abstract type AbstractTermLike{T} <: AbstractPolynomialLike{T} end
26+
2327
"""
2428
AbstractMonomialLike
2529
@@ -55,6 +59,7 @@ abstract type AbstractTerm{T} <: AbstractTermLike{T} end
5559
Abstract type for a polynomial of coefficient type `T`, i.e. a sum of `AbstractTerm{T}`s.
5660
"""
5761
abstract type AbstractPolynomial{T} <: AbstractPolynomialLike{T} end
62+
MA.mutability(::Type{<:AbstractPolynomial}) = MA.IsMutable()
5863

5964
const APL{T} = AbstractPolynomialLike{T}
6065

‎src/conversion.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,13 +2,14 @@ export variable
22

33
function convertconstant end
44
Base.convert(::Type{P}, α) where P<:APL = convertconstant(P, α)
5-
function Base.convert(::Type{P}, p::P) where P<:AbstractPolynomial
5+
function Base.convert(::Type{P}, p::P) where {T, P<:AbstractPolynomial{T}}
66
return p
77
end
88
function Base.convert(::Type{P}, p::APL) where {T, P<:AbstractPolynomial{T}}
99
return convert(P, polynomial(p, T))
1010
end
1111

12+
MA.scaling(p::AbstractPolynomialLike{T}) where {T} = convert(T, p)
1213
Base.convert(::Type{Any}, p::APL) = p
1314
# Conversion polynomial -> scalar
1415
function Base.convert(S::Type{<:Union{Number, T}}, p::APL{T}) where T

‎src/operators.jl‎

Lines changed: 196 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -20,15 +20,147 @@ end
2020
for op in [:+, :-, :(==)]
2121
@eval Base.$op(p1::APL, p2::APL) = $op(promote(p1, p2)...)
2222
end
23+
# Promotion between `I` and `1` is `Any`.
24+
# Promotion between `I` and `2I` is `UniformScaling`.
25+
for op in [:+, :-]
26+
@eval Base.$op(p1::APL, p2::APL{<:LinearAlgebra.UniformScaling}) = $op(p1, mapcoefficientsnz(J -> J.λ, p2))
27+
@eval Base.$op(p1::APL{<:LinearAlgebra.UniformScaling}, p2::APL) = $op(mapcoefficientsnz(J -> J.λ, p1), p2)
28+
@eval Base.$op(p1::APL{<:LinearAlgebra.UniformScaling}, p2::APL{<:LinearAlgebra.UniformScaling}) = $op(mapcoefficientsnz(J -> J.λ, p1), p2)
29+
end
2330
Base.isapprox(p1::APL, p2::APL; kwargs...) = isapprox(promote(p1, p2)...; kwargs...)
2431

2532
# @eval $op(p::APL, α) = $op(promote(p, α)...) would be less efficient
2633
for (op, fun) in [(:+, :plusconstant), (:-, :minusconstant), (:*, :multconstant), (:(==), :eqconstant)]
2734
@eval Base.$op(p::APL, α) = $fun(p, α)
2835
@eval Base.$op(α, p::APL) = $fun(α, p)
2936
end
37+
# Fix ambiguity between above methods and methods in MA
38+
Base.:+(::MA.Zero, p::APL) = MA.copy_if_mutable(p)
39+
Base.:+(p::APL, ::MA.Zero) = MA.copy_if_mutable(p)
40+
constant_function(::typeof(+)) = plusconstant
41+
constant_function(::typeof(-)) = minusconstant
42+
MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::APL, α) = MA.mutable_operate!(constant_function(op), p, α)
43+
MA.mutable_operate_to!(output::AbstractPolynomial, op::typeof(*), α, p::APL) = MA.mutable_operate_to!(output, multconstant, α, p)
44+
MA.mutable_operate_to!(output::AbstractPolynomial, op::typeof(*), p::APL, α) = MA.mutable_operate_to!(output, multconstant, p, α)
45+
46+
function polynomial_merge!(
47+
n1::Int, n2::Int, get1::Function, get2::Function,
48+
set::Function, push::Function, compare_monomials::Function,
49+
combine::Function, keep::Function, resize::Function)
50+
buffer = nothing
51+
i = j = k = 1
52+
# Invariant:
53+
# The terms p[0] -> p[k-1] are sorted and are smaller than the remaining terms.
54+
# The terms p[k] -> p[i-1] are garbage.
55+
# The terms p[i] -> p[end] are sorted and still need to be added.
56+
# The terms q[j] -> p[end] are sorted and still need to be added.
57+
# If `buffer` is not empty:
58+
# The terms in `buffer` are sorted and still need to be added.
59+
# Moreover, they are smaller than the terms p[i] -> p[end].
60+
while i <= n1 && j <= n2
61+
@assert buffer === nothing || isempty(buffer)
62+
comp = compare_monomials(i, j)
63+
if comp > 0
64+
if k == i
65+
t0 = get1(i)
66+
if buffer === nothing
67+
buffer = DataStructures.Queue{typeof(t0)}()
68+
end
69+
DataStructures.enqueue!(buffer, t0)
70+
i += 1
71+
end
72+
set(k, get2(j))
73+
j += 1
74+
k += 1
75+
elseif iszero(comp)
76+
combine(i, j)
77+
if keep(i)
78+
if k != i
79+
@assert k < i
80+
set(k, get1(i))
81+
end
82+
k += 1
83+
end
84+
i += 1
85+
j += 1
86+
else
87+
if k != i
88+
set(k, get1(i))
89+
end
90+
i += 1
91+
k += 1
92+
end
93+
while buffer !== nothing && !isempty(buffer) && j <= n2
94+
@assert i == k
95+
t = DataStructures.front(buffer)
96+
comp = compare_monomials(t, j)
97+
if comp >= 0
98+
if comp > 0
99+
t = get2(j)
100+
else
101+
t = combine(t, j)
102+
end
103+
j += 1
104+
end
105+
if comp <= 0
106+
DataStructures.dequeue!(buffer)
107+
end
108+
# if `comp` is zero, we called `combine` so `t`
109+
# might not be kept. If `comp` is not zero, we
110+
# skip the `keep` call that might be costly.
111+
if iszero(comp) && !keep(t)
112+
continue
113+
end
114+
if k <= n1
115+
DataStructures.enqueue!(buffer, get1(i))
116+
set(k, t)
117+
else
118+
push(t)
119+
n1 += 1
120+
end
121+
i += 1
122+
k += 1
123+
end
124+
end
125+
if buffer !== nothing && !isempty(buffer)
126+
@assert j == n2 + 1
127+
@assert i == k
128+
n = length(buffer)
129+
resize(n1 + n)
130+
for k in n1:-1:i
131+
set(k + n, get1(k))
132+
end
133+
for k in i:(i + n - 1)
134+
set(k, DataStructures.dequeue!(buffer))
135+
end
136+
n1 += n
137+
else
138+
len = (k - 1) + (n2 - (j - 1)) + (n1 - (i - 1))
139+
if n1 < len
140+
resize(len)
141+
end
142+
while j <= n2
143+
set(k, get2(j))
144+
j += 1
145+
k += 1
146+
end
147+
@assert j == n2 + 1
148+
while i <= n1
149+
set(k, get1(i))
150+
i += 1
151+
k += 1
152+
end
153+
if len < n1
154+
resize(len)
155+
end
156+
@assert i == n1 + 1
157+
@assert k == len + 1
158+
end
159+
return
160+
end
161+
30162

31-
MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomial) = MA.mutable_operate_to!(p, op, p, q)
163+
#MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomial) = MA.mutable_operate_to!(p, op, p, q)
32164
MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomialLike) = MA.mutable_operate!(op, p, polynomial(q))
33165

34166
# Special case AbstractArrays of APLs
@@ -41,6 +173,18 @@ for op in [:+, :-, :*]
41173
end
42174
Base.:/(A::AbstractArray{<:APL}, p::APL) = map(f -> f / p, A)
43175

176+
function mul_to_terms!(ts::Vector{<:AbstractTerm}, p1::APL, p2::APL)
177+
for t1 in terms(p1)
178+
for t2 in terms(p2)
179+
push!(ts, t1 * t2)
180+
end
181+
end
182+
return ts
183+
end
184+
function Base.:*(p::AbstractPolynomial, q::AbstractPolynomial)
185+
polynomial(mul_to_terms!(MA.promote_operation(*, termtype(p), termtype(q))[], p, q))
186+
end
187+
44188
Base.isapprox(p::APL, α; kwargs...) = isapprox(promote(p, α)...; kwargs...)
45189
Base.isapprox(α, p::APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)
46190

@@ -52,10 +196,16 @@ Base.:*(p::Union{APL, RationalPoly}) = p
52196

53197
# Avoid adding a zero constant that might artificially increase the Newton polytope
54198
# Need to add polynomial conversion for type stability
55-
plusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, Base.promote_op(+, S, T)) : p + constantterm(α, p)
56-
plusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial( p, Base.promote_op(+, S, T)) : constantterm(α, p) + p
57-
minusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, Base.promote_op(-, S, T)) : p - constantterm(α, p)
58-
minusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial(-p, Base.promote_op(-, S, T)) : constantterm(α, p) - p
199+
plusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(+, S, T)) : p + constantterm(α, p)
200+
plusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(+, S, T)) : constantterm(α, p) + p
201+
function MA.mutable_operate!(::typeof(plusconstant), p::APL, α)
202+
if !iszero(α)
203+
MA.mutable_operate!(+, p, constantterm(α, p))
204+
end
205+
return p
206+
end
207+
minusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(-, S, T)) : p - constantterm(α, p)
208+
minusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial(-p, MA.promote_operation(-, S, T)) : constantterm(α, p) - p
59209

60210
# Coefficients and variables commute
61211
multconstant(α, v::AbstractVariable) = multconstant(α, monomial(v)) # TODO linear term
@@ -64,7 +214,7 @@ multconstant(m::AbstractMonomialLike, α) = multconstant(α, m)
64214
_multconstant(α, f, t::AbstractTermLike) = mapcoefficientsnz(f, t)
65215
function _multconstant(α::T, f, p::AbstractPolynomial{S}) where {S, T}
66216
if iszero(α)
67-
zero(polynomialtype(p, Base.promote_op(*, T, S)))
217+
zero(polynomialtype(p, MA.promote_operation(*, T, S)))
68218
else
69219
mapcoefficientsnz(f, p)
70220
end
@@ -74,6 +224,22 @@ _multconstant(α, f, p::AbstractPolynomialLike) = _multconstant(α, f, polynomia
74224
multconstant(α, p::AbstractPolynomialLike) = _multconstant(α, β -> α*β, p)
75225
multconstant(p::AbstractPolynomialLike, α) = _multconstant(α, β -> β*α, p)
76226

227+
function mapcoefficientsnz_to! end
228+
229+
function _multconstant_to!(output, α, f, p)
230+
if iszero(α)
231+
MA.mutable_operate!(zero, output)
232+
else
233+
mapcoefficientsnz_to!(output, f, p)
234+
end
235+
end
236+
function MA.mutable_operate_to!(output, ::typeof(multconstant), p::APL, α)
237+
_multconstant_to!(output, α, β -> β*α, p)
238+
end
239+
function MA.mutable_operate_to!(output, ::typeof(multconstant), α, p::APL)
240+
_multconstant_to!(output, α, β -> α*β, p)
241+
end
242+
77243
MA.mutable_operate_to!(output::AbstractMonomial, ::typeof(*), m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents_to!(output, +, m1, m2)
78244
MA.mutable_operate!(::typeof(*), m1::AbstractMonomial, m2::AbstractMonomialLike) = mapexponents!(+, m1, m2)
79245
Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(+, m1, m2)
@@ -100,9 +266,9 @@ end
100266
for op in [:+, :-]
101267
@eval begin
102268
Base.$op(t1::AbstractTermLike, t2::AbstractTermLike) = $op(term(t1), term(t2))
103-
Base.$op(t1::AbstractTerm, t2::AbstractTerm) = $op(promote(t1, t2)...)
269+
Base.$op(t1::AbstractTerm, t2::AbstractTerm) = $op(_promote_terms(t1, t2)...)
104270
function Base.$op(t1::TT, t2::TT) where {T, TT <: AbstractTerm{T}}
105-
S = Base.promote_op($op, T, T)
271+
S = MA.promote_operation($op, T, T)
106272
# t1 > t2 would compare the coefficient in case the monomials are equal
107273
# and it will throw a MethodError in case the coefficients are not comparable
108274
if monomial(t1) == monomial(t2)
@@ -115,6 +281,18 @@ for op in [:+, :-]
115281
end
116282
end
117283
end
284+
_promote_terms(t1, t2) = promote(t1, t2)
285+
# Promotion between `I` and `1` is `Any`.
286+
_promote_terms(t1::AbstractTerm, t2::AbstractTerm{<:LinearAlgebra.UniformScaling}) = _promote_terms(t1, coefficient(t2).λ * monomial(t2))
287+
_promote_terms(t1::AbstractTerm{<:LinearAlgebra.UniformScaling}, t2::AbstractTerm) = _promote_terms(coefficient(t1).λ * monomial(t1), t2)
288+
# Promotion between `I` and `2I` is `UniformScaling`, not `UniformScaling{Int}`.
289+
function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{S}}, t2::AbstractTerm{LinearAlgebra.UniformScaling{T}}) where {S<:Number, T<:Number}
290+
U = LinearAlgebra.UniformScaling{promote_type(S, T)}
291+
_promote_terms(MA.scaling_convert(U, coefficient(t1)) * monomial(t1), MA.scaling_convert(U, coefficient(t2)) * monomial(t2))
292+
end
293+
function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{T}}, t2::AbstractTerm{LinearAlgebra.UniformScaling{T}}) where T<:Number
294+
return promote(t1, t2)
295+
end
118296

119297
LinearAlgebra.adjoint(v::AbstractVariable) = v
120298
LinearAlgebra.adjoint(m::AbstractMonomial) = m
@@ -132,6 +310,9 @@ LinearAlgebra.dot(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike) = p1'
132310
LinearAlgebra.dot(x, p::AbstractPolynomialLike) = x' * p
133311
LinearAlgebra.dot(p::AbstractPolynomialLike, x) = p' * x
134312

313+
LinearAlgebra.symmetric_type(PT::Type{<:APL}) = PT
314+
LinearAlgebra.symmetric(p::APL, ::Symbol) = p
315+
135316
# Amazingly, this works! Thanks, StaticArrays.jl!
136317
"""
137318
Convert a tuple of variables into a static vector to allow array-like usage.
@@ -146,15 +327,15 @@ Base.:^(x::AbstractPolynomialLike, p::Integer) = Base.power_by_squaring(x, p)
146327
function MA.mutable_operate_to!(output::AbstractPolynomial, ::typeof(MA.add_mul), x, args::Vararg{Any, N}) where N
147328
return MA.mutable_operate_to!(output, +, x, *(args...))
148329
end
149-
function MA.mutable_operate!(::typeof(MA.add_mul), x, args::Vararg{Any, N}) where N
150-
return MA.mutable_operate!(+, x, *(args...))
330+
function MA.mutable_operate!(::typeof(MA.add_mul), x, y, z, args::Vararg{Any, N}) where N
331+
return MA.mutable_operate!(+, x, *(y, z, args...))
151332
end
152-
MA.buffer_for(::typeof(MA.add_mul), ::Type{<:AbstractPolynomial}, args::Vararg{Any, N}) where {N} = MA.promote_operation(*, args...)
153-
function MA.mutable_buffered_operate_to!(buffer::AbstractPolynomial, output::AbstractPolynomial, ::typeof(MA.add_mul), x, args::Vararg{Any, N}) where N
154-
product = MA.operate_to!(buffer, *, args...)
333+
MA.buffer_for(::typeof(MA.add_mul), ::Type{<:AbstractPolynomial}, args::Vararg{Type, N}) where {N} = zero(MA.promote_operation(*, args...))
334+
function MA.mutable_buffered_operate_to!(buffer::AbstractPolynomial, output::AbstractPolynomial, ::typeof(MA.add_mul), x, y, z, args::Vararg{Any, N}) where N
335+
product = MA.operate_to!(buffer, *, y, z, args...)
155336
return MA.mutable_operate_to!(output, +, x, product)
156337
end
157-
function MA.mutable_buffered_operate!(buffer::AbstractPolynomial, ::typeof(MA.add_mul), x, args::Vararg{Any, N}) where N
158-
product = MA.operate_to!(buffer, *, args...)
338+
function MA.mutable_buffered_operate!(buffer::AbstractPolynomial, ::typeof(MA.add_mul), x, y, z, args::Vararg{Any, N}) where N
339+
product = MA.operate_to!(buffer, *, y, z, args...)
159340
return MA.mutable_operate!(+, x, product)
160341
end

‎src/polynomial.jl‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,28 @@ function uniqterms(ts::AbstractVector{T}) where T <: AbstractTerm
104104
end
105105
result
106106
end
107+
function uniqterms!(ts::AbstractVector{<: AbstractTerm})
108+
i = firstindex(ts)
109+
for j in Iterators.drop(eachindex(ts), 1)
110+
if !iszero(ts[j])
111+
if monomial(ts[i]) == monomial(ts[j])
112+
ts[i] = MA.add!(coefficient(ts[i]), coefficient(ts[j])) * monomial(ts[i])
113+
else
114+
if !iszero(ts[i])
115+
i += 1
116+
end
117+
ts[i] = MA.copy_if_mutable(ts[j])
118+
end
119+
end
120+
end
121+
if i < length(ts)
122+
if iszero(ts[i])
123+
i -= 1
124+
end
125+
resize!(ts, i)
126+
end
127+
ts
128+
end
107129
polynomial(ts::AbstractVector{<:AbstractTerm}, s::SortedState) = polynomial(uniqterms(ts), SortedUniqState())
108130
polynomial(ts::AbstractVector{<:AbstractTerm}, s::UnsortedState=MessyState()) = polynomial(sort(ts, lt=(>)), sortstate(s))
109131

@@ -374,4 +396,5 @@ function Base.round(p::AbstractPolynomialLike; args...)
374396
polynomial(round.(terms(p); args...), SortedState())
375397
end
376398

399+
Base.ndims(::Type{<:AbstractPolynomialLike}) = 0
377400
Base.broadcastable(p::AbstractPolynomialLike) = Ref(p)

0 commit comments

Comments
 (0)