Skip to content

Commit 538c811

Browse files
blegatblegat
authored
Faster monomial multiplication (#87)
* Faster monomial multiplication * Refactor mutable one as well
1 parent f04c8b0 commit 538c811

File tree

3 files changed

+118
-98
lines changed

3 files changed

+118
-98
lines changed

‎src/cmult.jl‎

Lines changed: 84 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -6,111 +6,124 @@ function Base.:(*)(x::PolyVar{true}, y::PolyVar{true})
66
Monomial{true}(x > y ? [x,y] : [y,x], [1,1])
77
end
88
end
9-
function multiplyvar(v::Vector{PolyVar{true}}, x::PolyVar{true})
10-
i = findfirst(w->w <= x, v)
9+
function multiplyvar(v::Vector{PolyVar{true}}, x::PolyVar{true}, z)
10+
i = findfirst(Base.Fix2(<=, x), v)
1111
if (i != nothing && i > 0) && v[i] == x
12-
multiplyexistingvar(v, x, i)
12+
copy(v), multiplyexistingvar(i, z)
1313
else
14-
insertvar(v, x, (i == nothing || i == 0) ? length(v)+1 : i)
14+
j = (i == nothing || i == 0) ? length(v)+1 : i
15+
insertvar(v, x, j), insertvar(z, x, j)
1516
end
1617
end
1718
function Base.:(*)(x::PolyVar{true}, y::Monomial{true})
18-
w, updatez = multiplyvar(y.vars, x)
19-
Monomial{true}(w, updatez(y.z))
19+
w, z = multiplyvar(y.vars, x, y.z)
20+
Monomial{true}(w, z)
2021
end
2122
Base.:(*)(y::MonomialVector{true}, x::PolyVar{true}) = x * y
2223
function Base.:(*)(x::PolyVar{true}, y::MonomialVector{true})
23-
w, updatez = multiplyvar(y.vars, x)
24-
MonomialVector{true}(w, updatez.(y.Z))
24+
w, Z = multiplyvar(y.vars, x, y.Z)
25+
MonomialVector{true}(w, Z)
2526
end
26-
function multdivmono!(output_variables::Vector{PolyVar{true}},
27-
v::Vector{PolyVar{true}}, x::Monomial{true}, op)
27+
function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2::Vector{Int}) where {F<:Function}
28+
@. output = op(z1, z2)
29+
return
30+
end
31+
function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2::Vector{Int}, n) where {F<:Function}
32+
resize!(output, n)
33+
@. output = op(z1, z2)
34+
return
35+
end
36+
function _operate_exponents_to!(output::Vector{Int}, op::F, z1::Vector{Int}, z2::Vector{Int}, n, maps) where {F<:Function}
37+
resize!(output, n)
38+
I = maps[1]; i = 1; lI = length(I)
39+
J = maps[2]; j = 1; lJ = length(J)
40+
while i <= lI || j <= lJ
41+
if i > lI || (j <= lJ && J[j] < I[i])
42+
output[J[j]] = op(0, z2[j])
43+
j += 1
44+
elseif j > lJ || (i <= lI && I[i] < J[j])
45+
output[I[i]] = op(z1[I[i]], 0)
46+
i += 1
47+
else
48+
@assert I[i] == J[j]
49+
output[I[i]] = op(z1[I[i]], z2[j])
50+
i += 1
51+
j += 1
52+
end
53+
end
54+
return
55+
end
56+
# Not used yet
57+
#function _operate_exponents!(op::F, Z::Vector{Vector{Int}}, z2::Vector{Int}, args::Vararg{Any,N}) where {F<:Function,N}
58+
# return Vector{Int}[_operate_exponents!(op, z, z2, args...) for z in Z]
59+
#end
60+
function multdivmono!(output, output_variables::Vector{PolyVar{true}},
61+
v::Vector{PolyVar{true}}, x::Monomial{true}, op, z)
2862
if v == x.vars
2963
if output_variables == v
30-
updatez! = (output, z) -> begin
31-
@. output = op(z, x.z)
32-
return
33-
end
64+
_operate_exponents_to!(output, op, z, x.z)
3465
else
3566
resize!(output_variables, length(v))
3667
copyto!(output_variables, v)
37-
updatez! = (output, z) -> begin
38-
resize!(output, length(output_variables))
39-
@. output = op(z, x.z)
40-
return
41-
end
68+
_operate_exponents_to!(output, op, z, x.z, length(output_variables))
4269
end
4370
else
4471
maps = mergevars_to!(output_variables, [v, x.vars])
45-
n = length(v)
46-
updatez! = (output, z) -> begin
47-
resize!(output, length(output_variables))
48-
I = maps[1]; i = 1; lI = length(I)
49-
J = maps[2]; j = 1; lJ = length(J)
50-
while i <= lI || j <= lJ
51-
if i > lI || (j <= lJ && J[j] < I[i])
52-
output[J[j]] = op(0, x.z[j])
53-
j += 1
54-
elseif j > lJ || (i <= lI && I[i] < J[j])
55-
output[I[i]] = op(z[I[i]], 0)
56-
i += 1
57-
else
58-
@assert I[i] == J[j]
59-
output[I[i]] = op(z[I[i]], x.z[j])
60-
i += 1
61-
j += 1
62-
end
63-
end
72+
_operate_exponents_to!(output, op, z, x.z, length(output_variables), maps)
73+
end
74+
return
75+
end
76+
function _operate_exponents(op::F, z1::Vector{Int}, z2::Vector{Int}) where {F<:Function}
77+
return op.(z1, z2)
78+
end
79+
function _operate_exponents(op::F, z1::Vector{Int}, z2::Vector{Int}, n, maps) where {F<:Function}
80+
newz = zeros(Int, n)
81+
I = maps[1]; i = 1; lI = length(I)
82+
J = maps[2]; j = 1; lJ = length(J)
83+
while i <= lI || j <= lJ
84+
if i > lI || (j <= lJ && J[j] < I[i])
85+
newz[J[j]] = op(0, z2[j])
86+
j += 1
87+
elseif j > lJ || (i <= lI && I[i] < J[j])
88+
newz[I[i]] = op(z1[i], 0)
89+
i += 1
90+
else
91+
@assert I[i] == J[j]
92+
newz[I[i]] = op(z1[i], z2[j])
93+
i += 1
94+
j += 1
6495
end
6596
end
66-
return updatez!
97+
return newz
6798
end
68-
function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op)
99+
function _operate_exponents(op::F, Z::Vector{Vector{Int}}, z2::Vector{Int}, args::Vararg{Any,N}) where {F<:Function,N}
100+
return Vector{Int}[_operate_exponents(op, z, z2, args...) for z in Z]
101+
end
102+
function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op, z)
69103
if v == x.vars
70104
w = copy(v)
71-
updatez = z -> op.(z, x.z)
105+
z_new = _operate_exponents(op, z, x.z)
72106
else
73107
w, maps = mergevars([v, x.vars])
74-
updatez = z -> begin
75-
newz = zeros(Int, length(w))
76-
I = maps[1]; i = 1; lI = length(I)
77-
J = maps[2]; j = 1; lJ = length(J)
78-
while i <= lI || j <= lJ
79-
if i > lI || (j <= lJ && J[j] < I[i])
80-
newz[J[j]] = op(0, x.z[j])
81-
j += 1
82-
elseif j > lJ || (i <= lI && I[i] < J[j])
83-
newz[I[i]] = op(z[i], 0)
84-
i += 1
85-
else
86-
@assert I[i] == J[j]
87-
newz[I[i]] = op(z[i], x.z[j])
88-
i += 1
89-
j += 1
90-
end
91-
end
92-
newz
93-
end
108+
z_new = _operate_exponents(op, z, x.z, length(w), maps)
94109
end
95-
w, updatez
110+
return w, z_new
96111
end
97112
function MP.mapexponents_to!(output::Monomial{true}, f::Function, x::Monomial{true}, y::Monomial{true})
98-
updatez! = multdivmono!(output.vars, x.vars, y, f)
99-
updatez!(output.z, x.z)
100-
return x
113+
multdivmono!(output.z, output.vars, x.vars, y, f, x.z)
114+
return output
101115
end
102116
function MP.mapexponents!(f::Function, x::Monomial{true}, y::Monomial{true})
103-
updatez! = multdivmono!(x.vars, x.vars, y, f)
104-
updatez!(x.z, x.z)
117+
multdivmono!(x.z, x.vars, x.vars, y, f, x.z)
105118
return x
106119
end
107120
function MP.mapexponents(f::Function, x::Monomial{true}, y::Monomial{true})
108-
w, updatez = multdivmono(x.vars, y, f)
109-
Monomial{true}(w, updatez(x.z))
121+
w, z = multdivmono(x.vars, y, f, x.z)
122+
return Monomial{true}(w, z)
110123
end
111124
function Base.:(*)(x::Monomial{true}, y::MonomialVector{true})
112-
w, updatez = multdivmono(y.vars, x, +)
113-
MonomialVector{true}(w, updatez.(y.Z))
125+
w, Z = multdivmono(y.vars, x, +, y.Z)
126+
return MonomialVector{true}(w, Z)
114127
end
115128
Base.:(*)(y::MonomialVector{true}, x::Monomial{true}) = x * y
116129
Base.:(*)(x::Monomial{true}, y::PolyVar{true}) = y * x

‎src/mult.jl‎

Lines changed: 24 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
1-
function multiplyexistingvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
2-
updatez = z -> begin
3-
newz = copy(z)
4-
newz[i] += 1
5-
newz
6-
end
7-
copy(v), updatez
1+
function multiplyexistingvar(i::Int, z::Vector{Int}) where {C}
2+
newz = copy(z)
3+
newz[i] += 1
4+
return newz
5+
end
6+
function multiplyexistingvar(i::Int, Z::Vector{Vector{Int}}) where {C}
7+
return Vector{Int}[multiplyexistingvar(i, z) for z in Z]
88
end
99
function insertvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
1010
n = length(v)
@@ -15,14 +15,21 @@ function insertvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
1515
w[I] = v[I]
1616
w[i] = x
1717
w[K] = v[J]
18-
updatez = z -> begin
19-
newz = Vector{Int}(undef, n+1)
20-
newz[I] = z[I]
21-
newz[i] = 1
22-
newz[K] = z[J]
23-
newz
24-
end
25-
w, updatez
18+
return w
19+
end
20+
function insertvar(z::Vector{Int}, x::PolyVar, i::Int)
21+
n = length(z)
22+
I = 1:i-1
23+
J = i:n
24+
K = J.+1
25+
newz = Vector{Int}(undef, n+1)
26+
newz[I] = z[I]
27+
newz[i] = 1
28+
newz[K] = z[J]
29+
return newz
30+
end
31+
function insertvar(Z::Vector{Vector{Int}}, x::PolyVar, i::Int)
32+
return Vector{Int}[insertvar(z, x, i) for z in Z]
2633
end
2734

2835
include("cmult.jl")
@@ -153,8 +160,8 @@ function MA.mutable_operate!(::typeof(*), p::Polynomial{C}, q::Polynomial{C}) wh
153160
end
154161

155162
# Overwrite this method for monomial-like terms because
156-
# otherwise it would check `iszero(α)` and in that case
157-
# dismiss of the variable of `p` by performing
163+
# otherwise it would check `iszero(α)` and in that case
164+
# dismiss of the variable of `p` by performing
158165
# `operate_to!(zero, output :: Polynomial )` which only
159166
# respects the variables that are stored already
160167
function MP._multconstant_to!(output::Polynomial, α, f, p :: DMonomialLike)

‎src/ncmult.jl‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ function multiplyvar(v::Vector{PolyVar{false}}, z::Vector{Int}, x::PolyVar{false
1212
i -= 1
1313
end
1414
if i > 0 && v[i] == x
15-
multiplyexistingvar(v, x, i)
15+
return copy(v), multiplyexistingvar(i, z)
1616
else
1717
# ---->
1818
# \ |\ |\
@@ -40,9 +40,9 @@ function multiplyvar(v::Vector{PolyVar{false}}, z::Vector{Int}, x::PolyVar{false
4040
end
4141

4242
if i <= length(v) && v[i] == x
43-
multiplyexistingvar(v, x, i)
43+
return copy(v), multiplyexistingvar(i, z)
4444
else
45-
insertvar(v, x, i)
45+
return insertvar(v, x, i), insertvar(z, x, i)
4646
end
4747
end
4848
end
@@ -52,7 +52,7 @@ function multiplyvar(x::PolyVar{false}, v::Vector{PolyVar{false}}, z::Vector{Int
5252
i += 1
5353
end
5454
if i <= length(v) && v[i] == x
55-
multiplyexistingvar(v, x, i)
55+
return copy(v), multiplyexistingvar(i, z)
5656
else
5757
# <----
5858
# \ |\ |\
@@ -79,19 +79,19 @@ function multiplyvar(x::PolyVar{false}, v::Vector{PolyVar{false}}, z::Vector{Int
7979
i -= 1
8080
end
8181
if i > 0 && v[i] == x
82-
multiplyexistingvar(v, x, i)
82+
return copy(v), multiplyexistingvar(i, z)
8383
else
84-
insertvar(v, x, i+1)
84+
return insertvar(v, x, i+1), insertvar(z, x, i+1)
8585
end
8686
end
8787
end
8888
function Base.:(*)(x::PolyVar{false}, y::Monomial{false})
89-
w, updatez = multiplyvar(x, y.vars, y.z)
90-
Monomial{false}(w, updatez(y.z))
89+
w, z = multiplyvar(x, y.vars, y.z)
90+
Monomial{false}(w, z)
9191
end
9292
function Base.:(*)(y::Monomial{false}, x::PolyVar{false})
93-
w, updatez = multiplyvar(y.vars, y.z, x)
94-
Monomial{false}(w, updatez(y.z))
93+
w, z = multiplyvar(y.vars, y.z, x)
94+
Monomial{false}(w, z)
9595
end
9696

9797
function Base.:(*)(x::Monomial{false}, y::Monomial{false})

0 commit comments

Comments
 (0)