Skip to content

Commit 5dc9654

Browse files
blegatblegat
authored
Update with MultivariatePolynomials#bl/term (#113)
* Update with MultivariatePolynomials#bl/term * Add missing termtype * Fix * Fix * Add promotion for subs * Reverse monomial order * Do the change for polynomial_merge! * pop! * pop! * Fixes suggested by VSCode * Preserve variables on term to poly conversion * Update tests with reordering * Fix monomials * Fix * Remove duplicate _vec * Use master branch of MP in ci * 1.0 -> 1.6 in CI
1 parent 8b448b8 commit 5dc9654

File tree

13 files changed

+87
-74
lines changed

13 files changed

+87
-74
lines changed

‎.github/workflows/ci.yml

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ jobs:
1717
- version: '1'
1818
os: ubuntu-latest
1919
arch: x64
20-
- version: '1.0'
20+
- version: '1.6'
2121
os: ubuntu-latest
2222
arch: x64
2323
- version: '1'
@@ -42,6 +42,11 @@ jobs:
4242
${{ runner.os }}-test-${{ env.cache-name }}-
4343
${{ runner.os }}-test-
4444
${{ runner.os }}-
45+
- name: MP#master
46+
shell: julia --project=@. {0}
47+
run: |
48+
using Pkg
49+
Pkg.add(PackageSpec(name="MultivariatePolynomials", rev="master"))
4550
- uses: julia-actions/julia-buildpkg@v1
4651
- uses: julia-actions/julia-runtest@v1
4752
with:

‎src/DynamicPolynomials.jl

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -32,10 +32,7 @@ MP.monomialtype(::PolyType{C}) where C = Monomial{C}
3232
#end
3333
MP.termtype(::Union{TermPoly{C, T}, Type{<:TermPoly{C, T}}}) where {C, T} = Term{C, T}
3434
MP.termtype(::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{T}) where {C, T} = Term{C, T}
35-
MP.polynomial(p::PolyType) = Polynomial(p)
36-
function MP.polynomial(p::PolyType{C}, ::Type{T}) where {C, T}
37-
return convert(Polynomial{C, T}, p)
38-
end
35+
MP.termtype(::Type{Polynomial{C}}) where {C} = Term{C}
3936
MP.polynomialtype(::Type{Term{C}}) where {C} = Polynomial{C}
4037
MP.polynomialtype(::Type{Term{C, T}}) where {T, C} = Polynomial{C, T}
4138
MP.polynomialtype(::Type{T}, ::Type{<:DMonomialLike{C}}) where {T, C} = Polynomial{C, T}

‎src/cmult.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,10 +8,10 @@ function Base.:(*)(x::PolyVar{true}, y::PolyVar{true})
88
end
99
function multiplyvar(v::Vector{PolyVar{true}}, x::PolyVar{true}, z)
1010
i = findfirst(Base.Fix2(<=, x), v)
11-
if (i != nothing && i > 0) && v[i] == x
11+
if (i !== nothing && i > 0) && v[i] == x
1212
copy(v), multiplyexistingvar(i, z)
1313
else
14-
j = (i == nothing || i == 0) ? length(v)+1 : i
14+
j = (i === nothing || i == 0) ? length(v)+1 : i
1515
insertvar(v, x, j), insertvar(z, x, j)
1616
end
1717
end

‎src/comp.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -155,12 +155,12 @@ function Base.isapprox(p::Polynomial{C, S}, q::Polynomial{C, T};
155155
ztol::Real=iszero(atol) ? Base.rtoldefault(S, T, 0) : atol) where {C, S, T}
156156
i = j = 1
157157
while i <= length(p.x) || j <= length(q.x)
158-
if i > length(p.x) || (j <= length(q.x) && q.x[j] > p.x[i])
158+
if i > length(p.x) || (j <= length(q.x) && q.x[j] < p.x[i])
159159
if !isapproxzero(q.a[j], ztol=ztol)
160160
return false
161161
end
162162
j += 1
163-
elseif j > length(q.x) || p.x[i] > q.x[j]
163+
elseif j > length(q.x) || p.x[i] < q.x[j]
164164
if !isapproxzero(p.a[i], ztol=ztol)
165165
return false
166166
end

‎src/diff.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
function MP.differentiate(m::Monomial{C}, x::PolyVar{C}) where C
22
i = findfirst(isequal(x), _vars(m))
3-
if (i == nothing || i == 0) || m.z[i] == 0
3+
if (i === nothing || i == 0) || m.z[i] == 0
44
zeroterm(m)
55
else
66
z = copy(m.z)

‎src/monovec.jl

Lines changed: 26 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ struct MonomialVector{C} <: AbstractVector{Monomial{C}}
88
function MonomialVector{C}(vars::Vector{PolyVar{C}}, Z::Vector{Vector{Int}}) where {C}
99
@assert !C || issorted(vars, rev=true)
1010
@assert all(z -> length(z) == length(vars), Z)
11-
@assert issorted(Z, rev=true, lt=grlex)
11+
@assert issorted(Z, lt=grlex)
1212
new(vars, Z)
1313
end
1414
end
@@ -52,6 +52,10 @@ function Base.deleteat!(x::MonomialVector, i)
5252
deleteat!(x.Z, i)
5353
return x
5454
end
55+
function Base.pop!(x::MonomialVector)
56+
pop!(x.Z)
57+
return x
58+
end
5559

5660
Base.firstindex(x::MonomialVector) = firstindex(x.Z)
5761
Base.lastindex(x::MonomialVector) = lastindex(x.Z)
@@ -84,26 +88,21 @@ MP.emptymonovec(::Type{<:DMonoVecElemNonConstant{C}}) where {C} = MonomialVector
8488

8589
function fillZfordeg!(Z, n, deg, ::Type{Val{true}}, filter::Function, ::Int)
8690
z = zeros(Int, n)
87-
z[1] = deg
91+
z[end] = deg
8892
while true
8993
if filter(z)
9094
push!(Z, z)
9195
z = copy(z)
9296
end
93-
if z[end] == deg
97+
if z[1] == deg
9498
break
9599
end
96-
sum = 1
97-
for j in (n-1):-1:1
98-
if z[j] != 0
99-
z[j] -= 1
100-
z[j+1] += sum
101-
break
102-
else
103-
sum += z[j+1]
104-
z[j+1] = 0
105-
end
106-
end
100+
i = findfirst(i -> !iszero(z[i]), n:-1:2)
101+
j = (n:-1:2)[i]
102+
p = z[j]
103+
z[j] = 0
104+
z[end] = p - 1
105+
z[j-1] += 1
107106
end
108107
end
109108
function fillZrec!(Z, z, i, n, deg, filter::Function)
@@ -121,17 +120,19 @@ function fillZrec!(Z, z, i, n, deg, filter::Function)
121120
end
122121
function fillZfordeg!(Z, n, deg, ::Type{Val{false}}, filter::Function, maxdeg::Int)
123122
z = zeros(Int, maxdeg * n - maxdeg + 1)
123+
start = length(Z) + 1
124124
fillZrec!(Z, z, 1, n, deg, filter)
125+
reverse!(view(Z, start:length(Z)))
125126
end
126127
# List exponents in decreasing Graded Lexicographic Order
127128
function getZfordegs(n, degs::AbstractVector{Int}, ::Type{Val{C}}, filter::Function) where C
128129
Z = Vector{Vector{Int}}()
129130
# For non-commutative, lower degree need to create a vector of exponent as large as for the highest degree
130131
maxdeg = isempty(degs) ? 0 : maximum(degs)
131-
for deg in sort(degs, rev=true)
132+
for deg in sort(degs)
132133
fillZfordeg!(Z, n, deg, Val{C}, filter, maxdeg)
133134
end
134-
@assert issorted(Z, rev=true, lt=grlex)
135+
@assert issorted(Z, lt=grlex)
135136
Z
136137
end
137138

@@ -191,7 +192,7 @@ end
191192
MP.sortmonovec(X::MonomialVector) = (1:length(X), X)
192193
function _sortmonovec(X::DMonoVec{C}) where {C}
193194
allvars, Z = buildZvarsvec(PolyVar{C}, X)
194-
σ = sortperm(Z, rev=true, lt=grlex)
195+
σ = sortperm(Z, lt=grlex)
195196
allvars, Z, σ
196197
end
197198
function _removedups!(Z::Vector{Vector{Int}}, σ::Vector{Int})
@@ -210,7 +211,7 @@ end
210211

211212
function MonomialVector{C}(X::DMonoVec{C}) where C
212213
allvars, Z = buildZvarsvec(PolyVar{C}, X)
213-
sort!(Z, rev=true, lt=grlex)
214+
sort!(Z, lt=grlex)
214215
dups = findall(i -> Z[i] == Z[i-1], 2:length(Z))
215216
deleteat!(Z, dups)
216217
MonomialVector{C}(allvars, Z)
@@ -231,21 +232,21 @@ function MP.mergemonovec(ms::Vector{MonomialVector{C}}) where {C}
231232
L = length.(ms)
232233
X = Vector{Monomial{C}}()
233234
while any(I .<= L)
234-
max = nothing
235+
min_monomial = nothing
235236
for i in 1:m
236237
if I[i] <= L[i]
237238
x = ms[i][I[i]]
238-
if max === nothing || max < x
239-
max = x
239+
if min_monomial === nothing || min_monomial > x
240+
min_monomial = x
240241
end
241242
end
242243
end
243-
@assert max !== nothing
244+
@assert min_monomial !== nothing
244245
# to ensure that max is no more a union
245-
max === nothing && return X
246-
push!(X, max)
246+
min_monomial === nothing && return X
247+
push!(X, min_monomial)
247248
for i in 1:m
248-
if I[i] <= L[i] && max == ms[i][I[i]]
249+
if I[i] <= L[i] && min_monomial == ms[i][I[i]]
249250
I[i] += 1
250251
end
251252
end

‎src/ncmult.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -96,11 +96,11 @@ end
9696

9797
function Base.:(*)(x::Monomial{false}, y::Monomial{false})
9898
i = findlast(z -> z > 0, x.z)
99-
if i == nothing || i == 0
99+
if i === nothing || i == 0
100100
return y
101101
end
102102
j = findfirst(z -> z > 0, y.z)
103-
if j == nothing || j == 0
103+
if j === nothing || j == 0
104104
return x
105105
end
106106
if x.vars[i] == y.vars[j]

‎src/operators.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,12 +21,12 @@ function _plusorminus_to!(a::Vector{U}, Z::Vector{Vector{Int}}, op::Function, p:
2121
i = j = 1
2222
while i <= nterms(p) || j <= nterms(q)
2323
z = zeros(Int, nvars)
24-
if j > nterms(q) || (i <= nterms(p) && _getindex(p, i).x > _getindex(q, j).x)
24+
if j > nterms(q) || (i <= nterms(p) && _getindex(p, i).x < _getindex(q, j).x)
2525
t = _getindex(p, i)
2626
z[maps[1]] = t.x.z
2727
α = MA.scaling_convert(U, MA.copy_if_mutable(t.α))
2828
i += 1
29-
elseif i > nterms(p) || _getindex(q, j).x > _getindex(p, i).x
29+
elseif i > nterms(p) || _getindex(q, j).x < _getindex(p, i).x
3030
t = _getindex(q, j)
3131
z[maps[2]] = t.x.z
3232
α = MA.scaling_convert(U, MA.operate(op, t.α))

‎src/poly.jl

Lines changed: 29 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -15,9 +15,20 @@ struct Polynomial{C, T} <: AbstractPolynomial{T}
1515
return p
1616
end
1717
end
18+
function Polynomial{C,T}(terms::AbstractVector{<:Term{C}}) where {C,T}
19+
a = T[coefficient(t) for t in terms]
20+
monos = Monomial{C}[monomial(t) for t in terms]
21+
allvars, Z = buildZvarsvec(PolyVar{C}, monos)
22+
x = MonomialVector{C}(allvars, Z)
23+
return Polynomial{C,T}(a, x)
24+
end
1825

1926
iscomm(::Type{Polynomial{C, T}}) where {C, T} = C
2027

28+
function _zero_with_variables(::Type{Polynomial{C, T}}, vars::Vector{PolyVar{C}}) where {C,T}
29+
return Polynomial(T[], MonomialVector{C}(vars, Vector{Int}[]))
30+
end
31+
2132
Base.broadcastable(p::Polynomial) = Ref(p)
2233
MA.mutable_copy(p::Polynomial{C, T}) where {C, T} = Polynomial{C, T}(MA.mutable_copy(p.a), MA.mutable_copy(p.x))
2334
Base.copy(p::Polynomial) = MA.mutable_copy(p)
@@ -35,29 +46,23 @@ Polynomial(af::Union{Function, Vector}, x::DMonoVec{C}) where {C} = Polynomial{C
3546
Polynomial{C, T}(p::Polynomial{C, T}) where {C, T} = p
3647

3748
Base.convert(::Type{Polynomial{C, T}}, p::Polynomial{C, T}) where {C, T} = p
38-
function Base.convert(::Type{Polynomial{C, T}},
39-
p::Polynomial{C, S}) where {C, S, T}
40-
return Polynomial{C}(convert(Vector{T}, p.a), p.x)
41-
end
42-
#function convert(::Type{Polynomial{C, T}},
43-
# p::AbstractPolynomialLike) where {C, T}
44-
# return convert(Polynomial{C, T}, polynomial(p, T))
45-
#end
46-
function Base.convert(::Type{Polynomial{C, T}}, t::Term{C}) where {C, T}
47-
return Polynomial{C, T}(T[t.α], [t.x])
48-
end
49-
function Base.convert(::Type{Polynomial{C, T}}, m::DMonomialLike{C}) where {C, T}
50-
return Polynomial(convert(Term{C, T}, m))
49+
function Base.convert(::Type{Polynomial{C, T}}, t::AbstractTermLike) where {C, T}
50+
if iszero(t)
51+
_zero_with_variables(Polynomial{C,T}, variables(t))
52+
else
53+
# `exponents(::PolyVar)` gives a tuple
54+
return Polynomial{C, T}([coefficient(t)], MonomialVector{C}(variables(t), [_vec(exponents(t))]))
55+
end
5156
end
52-
function MP.convertconstant(::Type{Polynomial{C, T}}, α) where {C, T}
53-
return Polynomial(convert(Term{C, T}, α))
57+
function Base.convert(::Type{Polynomial{C, T}}, p::AbstractPolynomialLike) where {C, T}
58+
return Polynomial{C, T}(terms(p))
5459
end
5560

5661
Polynomial{C}(p::Union{Polynomial{C}, Term{C}, Monomial{C}, PolyVar{C}}) where {C} = Polynomial(p)
5762
Polynomial{C}(α) where {C} = Polynomial(Term{C}(α))
5863

5964
Polynomial(p::Polynomial) = p
60-
Polynomial(t::Term{C, T}) where {C, T} = Polynomial{C, T}([t.α], [t.x])
65+
Polynomial(t::Term{C, T}) where {C, T} = convert(Polynomial{C, T}, mutable_copy(t))
6166
Polynomial(x::Union{PolyVar{C}, Monomial{C}}) where {C} = Polynomial(Term{C}(x))
6267

6368
#Base.convert(::Type{TermContainer{C, T}}, p::Polynomial{C}) where {C, T} = Polynomial{C, T}(p)
@@ -123,24 +128,24 @@ MP.extdegree(p::Polynomial) = extdegree(p.x)
123128
MP.mindegree(p::Polynomial) = mindegree(p.x)
124129
MP.maxdegree(p::Polynomial) = maxdegree(p.x)
125130

126-
MP.leadingcoefficient(p::Polynomial{C, T}) where {C, T} = iszero(p) ? zero(T) : first(p.a)
127-
MP.leadingmonomial(p::Polynomial) = iszero(p) ? constantmonomial(p) : first(p.x)
128-
MP.leadingterm(p::Polynomial) = iszero(p) ? zeroterm(p) : first(terms(p))
131+
MP.leadingcoefficient(p::Polynomial{C, T}) where {C, T} = iszero(p) ? zero(T) : last(p.a)
132+
MP.leadingmonomial(p::Polynomial) = iszero(p) ? constantmonomial(p) : last(p.x)
133+
MP.leadingterm(p::Polynomial) = iszero(p) ? zeroterm(p) : last(terms(p))
129134

130135
function MP.removeleadingterm(p::Polynomial)
131-
Polynomial(p.a[2:end], p.x[2:end])
136+
Polynomial(p.a[1:end-1], p.x[1:end-1])
132137
end
133138
function MA.operate!(::typeof(MP.removeleadingterm), p::Polynomial)
134-
deleteat!(p.a, 1)
135-
deleteat!(p.x, 1)
139+
pop!(p.a)
140+
pop!(p.x)
136141
return p
137142
end
138143
function MP.removemonomials(p::Polynomial, x::MonomialVector)
139144
# use the fact that monomials are sorted to do this O(n) instead of O(n^2)
140145
j = 1
141146
I = Int[]
142147
for (i,t) in enumerate(p)
143-
while j <= length(x) && x[j] > t.x
148+
while j <= length(x) && x[j] < t.x
144149
j += 1
145150
end
146151
if j > length(x) || x[j] != t.x
@@ -166,7 +171,7 @@ end
166171

167172
function removedups_to!(a::Vector{T}, Z::Vector{Vector{Int}},
168173
adup::Vector{T}, Zdup::Vector{Vector{Int}}) where T
169-
σ = sortperm(Zdup, rev=true, lt=grlex)
174+
σ = sortperm(Zdup, lt=grlex)
170175
i = 0
171176
j = 1
172177
while j <= length(adup)

‎src/subs.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,6 +101,11 @@ function _subs(::MP.Subs, p::Polynomial{C, T}, vals::AbstractVector{S}) where {C
101101
return q
102102
end
103103

104+
function MA.promote_operation(::typeof(MP.substitute), ::Type{MP.Subs}, ::Type{Monomial{C}}, ::Type{Pair{PolyVar{C},T}}) where {C,T}
105+
U = MA.promote_operation(^, T, Int)
106+
return Term{C,U}
107+
end
108+
104109
function MP.substitute(st::MP.AbstractSubstitutionType, p::PolyType, s::MP.Substitutions)
105110
_subs(st, p, subsmap(st, _vars(p), s))
106111
end

0 commit comments

Comments
 (0)