Merge pull request #52 from JuliaAlgebra/bl/mutable_arithmetics

Implement MutableArithmetics

Author: blegat

Project.toml

@@ -4,13 +4,17 @@ repo = "https://github.com/JuliaAlgebra/DynamicPolynomials.jl.git"
 version = "0.3.3"
 
 [deps]
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Future = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-MultivariatePolynomials = "~0.3.2"
+MultivariatePolynomials = "~0.3.3"
+MutableArithmetics = "0.1.1"
 Reexport = ">= 0.2"
 julia = "1"

src/DynamicPolynomials.jl

@@ -1,16 +1,25 @@
 module DynamicPolynomials
 
+import Future # For `copy!`
+
 using Reexport
 @reexport using MultivariatePolynomials
 const MP = MultivariatePolynomials
 
+import MutableArithmetics
+const MA = MutableArithmetics
+
+using DataStructures
 
 include("var.jl")
 include("mono.jl")
 const DMonomialLike{C} = Union{Monomial{C}, PolyVar{C}}
+MA.mutability(::Type{<:Monomial}) = MA.IsMutable()
 include("term.jl")
+MA.mutability(::Type{Term{C, T}}) where {C, T} = MA.mutability(T)
 include("monovec.jl")
 include("poly.jl")
+MA.mutability(::Type{<:Polynomial}) = MA.IsMutable()
 const TermPoly{C, T} = Union{Term{C, T}, Polynomial{C, T}}
 const PolyType{C} = Union{Polynomial{C}, Term{C}, Monomial{C}, PolyVar{C}}
 MP.variable_union_type(::Union{PolyType{C}, Type{<:PolyType{C}}}) where {C} = PolyVar{C}
@@ -35,6 +44,7 @@ MP.variables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = _v
 MP.nvariables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = length(_vars(p))
 MP.similarvariable(p::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{Val{V}}) where {C, V} = PolyVar{C}(string(V))
 MP.similarvariable(p::Union{PolyType{C}, Type{<:PolyType{C}}}, V::Symbol) where {C} = PolyVar{C}(string(V))
+
 include("promote.jl")
 
 include("operators.jl")

src/cmult.jl

@@ -23,6 +23,38 @@ function Base.:(*)(x::PolyVar{true}, y::MonomialVector{true})
     w, updatez = multiplyvar(y.vars, x)
     MonomialVector{true}(w, updatez.(y.Z))
 end
+function multdivmono!(output_variables::Vector{PolyVar{true}},
+                      v::Vector{PolyVar{true}}, x::Monomial{true}, op)
+    if v == x.vars
+        updatez! = (output, z) -> @. output = op(z, x.z)
+    else
+        w, maps = mergevars([v, x.vars])
+        n = length(v)
+        resize!(output_variables, length(w))
+        output_variables[maps[1]] = v[1:n]
+        updatez! = (output, z) -> begin
+            resize!(output, length(w))
+            z[maps[1]] = z[1:n]
+            I = maps[1]; i = 1; lI = length(I)
+            J = maps[2]; j = 1; lJ = length(J)
+            while i <= lI || j <= lJ
+                if i > lI || (j <= lJ && J[j] < I[i])
+                    output[J[j]] = op(0, x.z[j])
+                    j += 1
+                elseif j > lJ || (i <= lI && I[i] < J[j])
+                    output[I[i]] = op(z[I[i]], 0)
+                    i += 1
+                else
+                    @assert I[i] == J[j]
+                    output[I[i]] = op(z[I[i]], x.z[j])
+                    i += 1
+                    j += 1
+                end
+            end
+        end
+    end
+    return updatez!
+end
 function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op)
     if v == x.vars
         # /!\ no copy done here for efficiency, do not mess up with vars
@@ -53,7 +85,17 @@ function multdivmono(v::Vector{PolyVar{true}}, x::Monomial{true}, op)
     end
     w, updatez
 end
-function MP.mapexponents(f, x::Monomial{true}, y::Monomial{true})
+function MP.mapexponents_to!(output::Monomial{true}, f::Function, x::Monomial{true}, y::Monomial{true})
+    updatez! = multdivmono!(output.vars, x.vars, y, f)
+    updatez!(output.z, x.z)
+    return x
+end
+function MP.mapexponents!(f::Function, x::Monomial{true}, y::Monomial{true})
+    updatez! = multdivmono!(x.vars, x.vars, y, f)
+    updatez!(x.z, x.z)
+    return x
+end
+function MP.mapexponents(f::Function, x::Monomial{true}, y::Monomial{true})
     w, updatez = multdivmono(x.vars, y, f)
     Monomial{true}(w, updatez(x.z))
 end

src/comp.jl

@@ -30,6 +30,21 @@ Base.isless(x::PolyVar{C}, y::PolyVar{C}) where C = isless(y.id, x.id)
 # Comparison of Monomial
 
 # graded lex ordering
+function samevars_grlex(x::Vector{Int}, y::Vector{Int})
+    @assert length(x) == length(y)
+    degx = sum(x)
+    degy = sum(y)
+    if degx != degy
+        degx - degy
+    else
+        for i in eachindex(x)
+            if x[i] != y[i]
+                return x[i] - y[i]
+            end
+        end
+        return 0
+    end
+end
 function mycomp(x::Monomial{C}, y::Monomial{C}) where C
     degx = degree(x)
     degy = degree(y)

src/mono.jl

@@ -35,7 +35,8 @@ end
 Monomial(α::Number) = convert(Monomial{true}, α)
 
 Base.broadcastable(m::Monomial) = Ref(m)
-Base.copy(m::M) where {M<:Monomial} = M(m.vars, copy(m.z))
+MA.mutable_copy(m::M) where {M<:Monomial} = M(copy(m.vars), copy(m.z))
+Base.copy(m::Monomial) = MA.mutable_copy(m)
 
 # Generate canonical reperesentation by removing variables that are not used
 function canonical(m::Monomial)

src/monovec.jl

@@ -35,8 +35,10 @@ function Base.hash(m::MonomialVector, u::UInt)
     end
 end
 
-# /!\ vars not copied, do not mess with vars
-Base.copy(m::MV) where {MV<:MonomialVector} = MV(m.vars, copy(m.Z))
+# vars copied as it may be modifed by `_add_variables!`
+# `mutable_copy` recursively copies the vector or vector of integers.
+MA.mutable_copy(m::MV) where {MV<:MonomialVector} = MV(copy(m.vars), MA.mutable_copy(m.Z))
+Base.copy(m::MonomialVector) = MA.mutable_copy(m)
 function Base.getindex(x::MV, I) where {MV<:MonomialVector}
     MV(x.vars, x.Z[sort(I)])
 end
@@ -239,3 +241,16 @@ function MP.mergemonovec(ms::Vector{MonomialVector{C}}) where {C}
     # There is no duplicate by construction
     return MonomialVector{C}(buildZvarsvec(PolyVar{C}, X)...)
 end
+
+function _add_variables!(monos::MonomialVector{C}, allvars::Vector{PolyVar{C}}, map) where C
+    Future.copy!(monos.vars, allvars)
+    if !isempty(monos.Z)
+        tmp = similar(first(monos.Z))
+        for z in monos.Z
+            Future.copy!(tmp, z)
+            resize!(z, length(allvars))
+            fill!(z, 0)
+            z[map] = tmp
+        end
+    end
+end

src/mult.jl

@@ -31,6 +31,20 @@ include("ncmult.jl")
 
 MP.multconstant(α, x::Monomial)   = Term(α, x)
 MP.mapcoefficientsnz(f::Function, p::Polynomial) = Polynomial(f.(p.a), p.x)
+function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, t::MP.AbstractTermLike)
+    MP.mapcoefficientsnz_to!(output, f, polynomial(t))
+end
+function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, p::Polynomial)
+    resize!(output.a, length(p.a))
+    @. output.a = f(p.a)
+    Future.copy!(output.x.vars, p.x.vars)
+    # TODO reuse the part of `Z` that is already in `output`.
+    resize!(output.x.Z, length(p.x.Z))
+    for i in eachindex(p.x.Z)
+        output.x.Z[i] = copy(p.x.Z[i])
+    end
+    return output
+end
 
 # I do not want to cast x to TermContainer because that would force the promotion of eltype(q) with Int
 function Base.:(*)(x::DMonomialLike, p::Polynomial)
@@ -48,7 +62,7 @@ function Base.:(*)(p::Polynomial, x::DMonomialLike)
 end
 
 function _term_poly_mult(t::Term{C, S}, p::Polynomial{C, T}, op::Function) where {C, S, T}
-    U = Base.promote_op(op, S, T)
+    U = MA.promote_operation(op, S, T)
     if iszero(t)
         zero(Polynomial{C, U})
     else
@@ -67,35 +81,64 @@ end
 Base.:(*)(p::Polynomial, t::Term) = _term_poly_mult(t, p, (α, β) -> β * α)
 Base.:(*)(t::Term, p::Polynomial) = _term_poly_mult(t, p, *)
 _sumprod(a, b) = a * b + a * b
-function Base.:(*)(p::Polynomial{C, S}, q::Polynomial{C, T}) where {C, S, T}
-    U = Base.promote_op(_sumprod, S, T)
-    if iszero(p) || iszero(q)
-        zero(Polynomial{C, U})
+function _mul(::Type{T}, p::AbstractPolynomialLike, q::AbstractPolynomialLike) where T
+    return _mul(T, polynomial(p), polynomial(q))
+end
+function _mul(::Type{T}, p::Polynomial{true}, q::Polynomial{true}) where T
+    samevars = _vars(p) == _vars(q)
+    if samevars
+        allvars = _vars(p)
     else
-        samevars = _vars(p) == _vars(q)
-        if samevars
-            allvars = _vars(p)
-        else
-            allvars, maps = mergevars([_vars(p), _vars(q)])
-        end
-        N = length(p)*length(q)
-        Z = Vector{Vector{Int}}(undef, N)
-        a = Vector{U}(undef, N)
-        i = 0
-        for u in p
-            for v in q
-                if samevars
-                    z = u.x.z + v.x.z
-                else
-                    z = zeros(Int, length(allvars))
-                    z[maps[1]] += u.x.z
-                    z[maps[2]] += v.x.z
-                end
-                i += 1
-                Z[i] = z
-                a[i] = u.α * v.α
+        allvars, maps = mergevars([_vars(p), _vars(q)])
+    end
+    N = length(p) * length(q)
+    Z = Vector{Vector{Int}}(undef, N)
+    a = Vector{T}(undef, N)
+    i = 0
+    for u in p
+        for v in q
+            if samevars
+                z = u.x.z + v.x.z
+            else
+                z = zeros(Int, length(allvars))
+                z[maps[1]] += u.x.z
+                z[maps[2]] += v.x.z
             end
+            i += 1
+            Z[i] = z
+            a[i] = u.α * v.α
         end
-        polynomialclean(allvars, a, Z)
     end
+    return allvars, a, Z
+end
+function Base.:(*)(p::Polynomial{true, S}, q::Polynomial{true, T}) where {S, T}
+    if iszero(p) || iszero(q)
+        zero(MA.promote_operation(*, typeof(p), typeof(q)))
+    else
+        polynomialclean(_mul(MA.promote_operation(MA.add_mul, S, T), p, q)...)
+    end
+end
+function MA.mutable_operate_to!(p::Polynomial{false, T}, ::typeof(*), q1::MP.AbstractPolynomialLike, q2::MP.AbstractPolynomialLike) where T
+    if iszero(q1) || iszero(q2)
+        MA.mutable_operate!(zero, p)
+    else
+        ts = Term{false, T}[]
+        MP.mul_to_terms!(ts, q1, q2)
+        # TODO do better than create tmp
+        tmp = polynomial!(ts)
+        Future.copy!(p.a, tmp.a)
+        Future.copy!(p.x.vars, tmp.x.vars)
+        Future.copy!(p.x.Z, tmp.x.Z)
+        return p
+    end
+end
+function MA.mutable_operate_to!(p::Polynomial{true, T}, ::typeof(*), q1::MP.AbstractPolynomialLike, q2::MP.AbstractPolynomialLike) where T
+    if iszero(q1) || iszero(q2)
+        MA.mutable_operate!(zero, p)
+    else
+        polynomialclean_to!(p, _mul(T, q1, q2)...)
+    end
+end
+function MA.mutable_operate!(::typeof(*), p::Polynomial{C}, q::Polynomial{C}) where C
+    return MA.mutable_operate_to!(p, *, p, q)
 end