MultivariatePolynomials integration

Manifest.toml

@@ -1,6 +1,187 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.42.0"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.11"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[MultivariatePolynomials]]
+deps = ["DataStructures", "LinearAlgebra", "MutableArithmetics"]
+path = "/Users/scheme/src/julia/MultivariatePolynomials"
+uuid = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
+version = "0.4.4"
+
+[[MutableArithmetics]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "ba8c0f8732a24facba709388c74ba99dcbfdda1e"
+uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+version = "1.0.0"
+
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+
+[[OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
 [[SIMD]]
 git-tree-sha1 = "39e3df417a0dd0c4e1f89891a281f82f5373ea3b"
 uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
 version = "3.4.0"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

Project.toml

@@ -4,9 +4,11 @@ authors = ["Yingbo Ma <[email protected]> and Chris Elrod <[email protected]>"]
 version = "0.1.0"
 
 [deps]
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 
 [compat]
+MultivariatePolynomials = "0.4"
 SIMD = "3"
 julia = "1"
 

src/SIMDPolynomials.jl

@@ -1,5 +1,8 @@
 module SIMDPolynomials
 
+using MultivariatePolynomials
+const Term = MultivariatePolynomials.Term
+const Polynomial = MultivariatePolynomials.Polynomial
 export Uninomial, Uniterm, Poly
 export PackedMonomial, Monomial, Term, MPoly
 export terms, coeffs, content, contprim
@@ -15,11 +18,11 @@ end
 
 debugmode() = false
 
-include("interface.jl")
+#include("interface.jl")
 include("utils.jl")
 include("monomial.jl")
 include("packedmonomial.jl")
-include("mpoly.jl")
-include("poly.jl")
+#include("mpoly.jl")
+#include("poly.jl")
 
 end

src/mpoly.jl

@@ -1,7 +1,9 @@
 # default polynomial type
-struct MPoly{T} <: AbstractPolynomial
+#=
+struct MPoly{T} <: AbstractPolynomial{T}
     terms::T
 end
+=#
 #MPoly() = MPoly(EMPTY_TERMS)
 #MPoly(x::AbstractTerm) = MPoly([x])
 #MPoly(x::M) where {M<:AbstractMonomial} = MPoly(Term(x))

src/packedmonomial.jl

@@ -49,8 +49,9 @@ struct PackedMonomial{L,E,K} <: AbstractMonomial
     PackedMonomial{L,E}(bits::NTuple{K,UInt64}) where {L,E,K} = new{L,new_E(Val(E)),K}(bits)
     PackedMonomial{L,E,K}(bits::NTuple{K,UInt64}) where {L,E,K} = new{L,new_E(Val(E)),K}(bits)
 end
+Base.copy(m::PackedMonomial) = m
 nvariables(x::PackedMonomial{L}) where L = L
-function exps(m::PackedMonomial{L,E,K}) where {L,E,K}
+function MultivariatePolynomials.exponents(m::PackedMonomial{L,E,K}) where {L,E,K}
     tup = m.bits
     k = 0
     vpu = var_per_UInt64(Val(E))
@@ -148,6 +149,9 @@ end
 
 # the first chunk is the total degree
 Base.isone(x::PackedMonomial) = iszero(first(x.bits))
+function MultivariatePolynomials.grlex(x::T, y::T) where {T<:PackedMonomial}
+    x.bits < y.bits ? -1 : (x == y ? 0 : 1)
+end
 Base.isless(x::T, y::T) where {T<:PackedMonomial} = x.bits < y.bits
 @noinline overflowed_error() = throw(Base.OverflowError("Try increasing E."))
 @inline function check_zero_mask(::Val{E}) where {E}
@@ -179,12 +183,39 @@ function Base.:*(x::T, y::T) where {L,E,T<:PackedMonomial{L,E}}
     o != zero(UInt64) && overflowed_error()
     return T(xys)
 end
-function Base.:/(x::T, y::T) where {L,E,T<:PackedMonomial{L,E}}
+
+function MultivariatePolynomials.divides(y::T, x::T) where {L,E,T<:PackedMonomial{L,E}}
+    xys = _fmap(-, x.bits, y.bits)
+    o = reduce_tup(|, _fmap(Base.Fix2(zero_bits, Val(E)), xys))
+
+    #o != zero(UInt64) && overflowed_error()
+    #return T(xys), o != zero(UInt64)
+    return o == zero(UInt64)
+end
+
+MultivariatePolynomials.constantmonomial(::M) where {M<:PackedMonomial} = constantmonomial(M)
+function MultivariatePolynomials.constantmonomial(M::Type{<:PackedMonomial})
+    M()
+end
+
+# TODO: make it fast
+function MultivariatePolynomials.mapexponents!(f::F, m1::M, m2::M) where {F<:Function, L, E, M<:PackedMonomial{L,E}}
+    MultivariatePolynomials.mapexponents(f, m1, m2)
+end
+function MultivariatePolynomials.mapexponents(f::F, m1::M, m2::M) where {F<:Function, L, E, M<:PackedMonomial{L,E}}
+    e1 = exponents(m1)
+    e2 = exponents(m2)
+    ne = map(f, e1, e2)
+    prod(((i, e),)->PackedMonomial{L,E}(i-1)^e, enumerate(ne))
+end
+
+function MultivariatePolynomials._div(x::T, y::T) where {L,E,T<:PackedMonomial{L,E}}
     xys = _fmap(-, x.bits, y.bits)
     o = reduce_tup(|, _fmap(Base.Fix2(zero_bits, Val(E)), xys))
 
     #o != zero(UInt64) && overflowed_error()
-    return T(xys), o != zero(UInt64)
+    #return T(xys), o != zero(UInt64)
+    return T(xys)#, o != zero(UInt64)
 end
 function Base.:^(x::T, y::Integer) where {L,E,T<:PackedMonomial{L,E}}
     xys = _fmap(*, x.bits, y)
@@ -263,3 +294,36 @@ function Base.show(io::IO, m::PackedMonomial{L,E}) where {L,E}
         end
     end
 end
+
+struct Variable{L,E} <: AbstractVariable
+    id::UInt
+end
+
+function Base.:(==)(v1::V, v2::V) where {V<:Variable}
+    v1 === v2
+end
+function Base.isless(v1::V, v2::V) where {V<:Variable}
+    isless(v1.id, v2.id)
+end
+function MultivariatePolynomials.name_base_indices(v::Variable)
+    "x", (v.id,)
+end
+MultivariatePolynomials.name(v::Variable) = string("x", int2lowerscript(v.id))
+MultivariatePolynomials.variable_union_type(::Type{<:PackedMonomial{L,E}}) where {L,E} = Variable{L,E}
+
+MultivariatePolynomials.variables(::Type{<:PackedMonomial{L,E}}) where {L, E} = ntuple(i->Variable{L,E}(i-1), Val(L))
+MultivariatePolynomials.variables(::M) where {M<:PackedMonomial} = variables(M)
+
+MultivariatePolynomials.termtype(M::Type{<:PackedMonomial}, ::Type{T}) where {T} = Term{T, M}
+MultivariatePolynomials.polynomialtype(::Type{Term{T, M}}) where {T, M<:PackedMonomial} = Polynomial{T, Term{T, M}, Vector{Term{T, M}}}
+
+MultivariatePolynomials.nvariables(p::Polynomial{T, Term{T, M}}) where {T, L, M<:PackedMonomial{L}} = L
+MultivariatePolynomials.variables(p::Polynomial{T, Term{T, M}}) where {T, M<:PackedMonomial} = variables(M)
+Base.:(^)(x::Variable{L,E}, p::Integer) where {L,E} = PackedMonomial{L,E}(x.id) ^ p
+function Base.:(==)(m1::T, m2::T) where {T<:PackedMonomial}
+    m1 === m2
+end
+
+function MultivariatePolynomials.substitute(::MultivariatePolynomials.Subs, v::V, p::Pair{V, Int64}) where {L,E,V<:Variable{L, E}}
+    v == p[1] ? p[2] : v
+end