diff --git a/.gitignore b/.gitignore index 46bfd1a..53e36cf 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ /docs/build/ statprof/ profile.pb.gz +Manifest.toml diff --git a/Manifest.toml b/Manifest.toml deleted file mode 100644 index 56acbde..0000000 --- a/Manifest.toml +++ /dev/null @@ -1,187 +0,0 @@ -# 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" diff --git a/Project.toml b/Project.toml index a5c8bac..98c8c11 100644 --- a/Project.toml +++ b/Project.toml @@ -5,10 +5,12 @@ version = "0.1.0" [deps] MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" SIMD = "fdea26ae-647d-5447-a871-4b548cad5224" [compat] -MultivariatePolynomials = "0.4" +MultivariatePolynomials = "0.5" +MutableArithmetics = "1" SIMD = "3" julia = "1" diff --git a/src/SIMDPolynomials.jl b/src/SIMDPolynomials.jl index 620b68a..28b9b8e 100644 --- a/src/SIMDPolynomials.jl +++ b/src/SIMDPolynomials.jl @@ -1,8 +1,9 @@ module SIMDPolynomials -using MultivariatePolynomials -const Term = MultivariatePolynomials.Term -const Polynomial = MultivariatePolynomials.Polynomial +import MutableArithmetics as MA +import MultivariatePolynomials as MP +const Term = MP.Term +const Polynomial = MP.Polynomial export Uninomial, Uniterm, Poly export PackedMonomial, Monomial, Term, MPoly export terms, coeffs, content, contprim @@ -18,11 +19,9 @@ end debugmode() = false -#include("interface.jl") +include("variable.jl") include("utils.jl") include("monomial.jl") include("packedmonomial.jl") -#include("mpoly.jl") -#include("poly.jl") end diff --git a/src/interface.jl b/src/interface.jl deleted file mode 100644 index 67e7be1..0000000 --- a/src/interface.jl +++ /dev/null @@ -1,165 +0,0 @@ -### -### Contract: mutation without explicit copy is undefined behavior. -### - -const PRETTY_PRINT = Ref(true) -const CoeffType = Union{Rational{<:Integer},Integer} - -### -### AbstractMonomial: isless, degree -### -abstract type AbstractMonomial <: Number end - -Base.isone(x::AbstractMonomial) = iszero(degree(x)) -Base.one(::Type{<:T}) where {T<:AbstractMonomial} = T() - -Base.:-(y::AbstractMonomial) = Term(-1, y) -Base.:+(y::T, x::T) where {T <: AbstractMonomial} = Term(y) + Term(x) -Base.:*(c::CoeffType, m::AbstractMonomial) = Term(c, m) -Base.:*(m::AbstractMonomial, c::CoeffType) = c * m -Base.:^(y::T, n::Integer) where {T <: AbstractMonomial} = iszero(n) ? T() : Base.power_by_squaring(y, n) -monomialtype(x::AbstractMonomial) = typeof(x) -nvariables(x::AbstractMonomial) = error("nvariables not implemented for $(typeof(x))") - -### -### AbstractTerm: `coeff` and `monomial` -### -abstract type AbstractTerm <: Number end -Base.promote_rule(::Type{C}, ::Type{M}) where {C<:CoeffType,M<:AbstractMonomial} = Term{C,M} -Base.promote_rule(t::Type{<:AbstractTerm}, ::Type{<:AbstractMonomial}) = t -Base.promote_rule(t::Type{<:AbstractTerm}, ::Type{<:CoeffType}) = t -emptyterm(T::Type{<:AbstractTerm}) = T[] - -degree(x::AbstractTerm) = degree(monomial(x)) -ismatch(x::T, y::T) where {T<:AbstractTerm} = monomial(x) == monomial(y) -Base.:(==)(x::AbstractTerm, y::AbstractTerm) = x === y || (coeff(x) == coeff(y) && monomial(x) == monomial(y)) -Base.copy(x::T) where {T<:AbstractTerm} = T(copy(coeff(x)), copy(monomial(x))) -Base.isless(x::T, y::T) where {T<:AbstractTerm} = isless(monomial(x), monomial(y)) -Base.iszero(x::AbstractTerm) = iszero(coeff(x)) -Base.isone(x::AbstractTerm) = isone(coeff(x)) && isone(monomial(x)) -Base.zero(t::T) where {T<:AbstractTerm} = parameterless_type(T)(zero(coeff(t)), one(monomial(t))) -Base.one(t::T) where {T<:AbstractTerm} = parameterless_type(T)(one(coeff(t)), one(monomial(t))) -Base.isinteger(x::AbstractTerm) = isinteger(coeff(x)) - -monomialtype(x::AbstractTerm) = monomialtype(monomial(x)) -exps(x::AbstractTerm) = exps(monomial(x)) - -Base.:*(x::T, y::T) where {T<:AbstractTerm} = T(coeff(x) * coeff(y), monomial(x) * monomial(y)) -# TODO -function Base.:(/)(y::T, x::T) where {T<:AbstractTerm} - m, fail = monomial(y) / monomial(x) - if fail - T(coeff(y), m), fail - else - T(coeff(y)/coeff(x), m), fail - end -end - -function Base.:+(x::T, y::T) where {T<:AbstractTerm} - if ismatch(x, y) - c = coeff(x) + coeff(y) - return iszero(c) ? MPoly(emptyterm(T)) : MPoly(T(c, monomial(x))) - else - if iszero(x) && iszero(y) - return MPoly(emptyterm(T)) - elseif iszero(x) - return MPoly(T[y]) - elseif iszero(y) - return MPoly(T[x]) - else - if x < y - x, y = y, x - end - return MPoly(T[x, y]) - end - end -end - -Base.:-(x::T) where {T<:AbstractTerm} = T(-coeff(x), monomial(x)) -function Base.:-(x::T, y::T) where {T<:AbstractTerm} - if ismatch(x, y) - c = coeff(x) - coeff(y) - return iszero(c) ? MPoly(T[]) : MPoly(T(c, monomial(x))) - else - y = -y - if x < y - x, y = y, x - end - return MPoly(T[x, y]) - end -end - -function Base.gcd(x::T, y::T) where {T<:AbstractTerm} - #g, a, b = gcd(monomial(x), monomial(y)) - g = gcd(monomial(x), monomial(y)) - gr = gcd(x.coeff, y.coeff) - return T(gr, g)#, Term(x.coeff / gr, a), Term(y.coeff / gr, b) -end - -print_coeff(io::IO, coeff) = isinteger(coeff) ? print(io, Integer(coeff)) : print(io, coeff) -function Base.show(io::IO, x::AbstractTerm) - printed = false - if coeff(x) != 1 - if coeff(x) == -1 - print(io, '-') - else - print_coeff(io, coeff(x)) - printed = true - PRETTY_PRINT[] || print(io, '*') - end - end - m = monomial(x) - if !(printed && isone(m)) - show(io, m) - end -end - -# default term type -struct Term{C,M<:AbstractMonomial} <: AbstractTerm - coeff::C - monomial::M -end -coeff(x::Term) = x.coeff -monomial(x::Term) = x.monomial -Term{M}(x) where {M<:AbstractMonomial} = Term(x, M()) -Term(x::M) where {M<:AbstractMonomial} = Term(1, x) -Term{A,B}(x::M) where {A,B,M<:AbstractMonomial} = Term(1, x) -Term{A,B}(x) where {A,B} = Term(x, B()) -#const EMPTY_TERM = Term[] -coefftype(::Type{<:Term{T}}) where T = T - -### -### AbstractPolynomial: terms, copy -### -abstract type AbstractPolynomial <: Number end - -coefftype(p::AbstractPolynomial) = coefftype(eltype(terms(p))) -Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:AbstractMonomial}) = p -Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:AbstractTerm}) = p -Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:CoeffType}) = p - -Base.iszero(x::AbstractPolynomial) = isempty(terms(x)) -Base.isone(x::AbstractPolynomial) = (ts = terms(x); length(ts) == 1 && isone(only(ts))) -Base.:(==)(x::T, y::T) where {T<:AbstractPolynomial} = x === y || (terms(x) == terms(y)) -monomialtype(p::AbstractPolynomial) = monomialtype(lt(p)) -monomials(p::AbstractPolynomial) = (monomial(t) for t in terms(p)) - -function Base.show(io::IO, p::AbstractPolynomial) - ts = terms(p) - if isempty(ts) - print(io, 0) - return - end - n = length(ts) - t1, tr = Iterators.peel(ts) - show(io, t1) - for t in tr - if t.coeff < 0 - print(io, " - ") - t = -t - else - print(io, " + ") - end - show(io, t) - end -end diff --git a/src/monomial.jl b/src/monomial.jl index c68f6be..9b2c6bd 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -2,10 +2,23 @@ const IDType = UInt32 const NOT_A_VAR = typemax(IDType) const EMPTY_IDS = IDType[] -struct Monomial <: AbstractMonomial +struct Variable <: AbstractVariable + id::IDType +end + +function Base.:(^)(x::Variable, p::Integer) + return Monomial([x.id for _ in 1:p]) +end + +function MP.monomial_type(::Type{Variable}) + return Monomial +end + +struct Monomial <: MP.AbstractMonomial ids::Vector{IDType} end Monomial() = Monomial(EMPTY_IDS) +MP.constant_monomial(::Union{Monomial,Type{Monomial}}) = Monomial() degree(x::Monomial) = length(x.ids) Base.copy(x::Monomial) = Monomial(copy(x.ids)) #TODO: optimize @@ -25,60 +38,13 @@ firstid(m::Monomial) = m.ids[1] degree(m::Monomial, id) = count(isequal(id), m.ids) rmid(m::Monomial, id) = Monomial(filter(!isequal(id), m.ids)) -const VARNAME_DICT = Dict{IDType, String}(0 => "x", 1 => "y", 2 => "z") -const SUPERSCRIPTS = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹'] -const LOWERSCRIPTS = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'] - -function int2superscript(x) - mapreduce(d->SUPERSCRIPTS[d + 1], *, Iterators.reverse(digits(x))) -end -function int2lowerscript(x) - mapreduce(d->LOWERSCRIPTS[d + 1], *, Iterators.reverse(digits(x))) -end -function print_single_monomial(io, v, o, star=false) - iszero(o) && return print(io, 1) - star && (PRETTY_PRINT[] || print(io, '*')) - print(io, VARNAME_DICT[v]) - if o > 1 - if PRETTY_PRINT[] - print(io, int2superscript(o)) - else - print(io, '^', o) - end - elseif o == 1 || o == 0 - else - error("unreachable") - end -end -function Base.show(io::IO, x::Monomial) - isempty(x.ids) && (print(io, '1'); return) - - ids = x.ids - v = first(ids) - count = 0 - star = false - for id in ids - if id == v - count += 1 - else - print_single_monomial(io, v, count, star) - star = true - v = id - count = 1 - end - end - if count > 0 - print_single_monomial(io, v, count, star) - end - return nothing -end - -function Base.:*(y::Monomial, x::Monomial) +function _p(y::Monomial, x::Monomial) ids = y.ids i = j = 1 n0 = length(ids) n1 = length(x.ids) - r = Monomial(similar(ids, n0 + n1)) + #r = Monomial(similar(ids, n0 + n1)) + r = Monomial(zeros(IDType, n0 + n1)) T = eltype(ids) M = typemax(T) for k in 1:(n0 + n1) @@ -95,30 +61,6 @@ function Base.:*(y::Monomial, x::Monomial) return r end -function Base.:(/)(y::Monomial, x::Monomial) - n = Monomial(similar(y.ids, 0)) - i = j = 1 - ids = y.ids - n0 = length(ids) - n1 = length(x.ids) - T = eltype(ids) - M = typemax(T) - while (i + j - 2) < (n0 + n1) - a = i <= n0 ? ids[i] : M - b = j <= n1 ? x.ids[j] : M - if a < b - push!(n.ids, a) - i += 1 - elseif a == b - i += 1 - j += 1 - else - return n, true - end - end - return n, false -end - Base.:(==)(x::Monomial, y::Monomial) = (x === y) || (x.ids == y.ids) # graded lex @@ -134,27 +76,116 @@ function Base.isless(x::Monomial, y::Monomial) end return false end +function MP.compare(x::Monomial, y::Monomial) + if x < y + return -1 + elseif x == y + return 0 + else + return 1 + end +end +MP.isconstant(x::Monomial) = isempty(x.ids) +MP.variables(x::Monomial) = Variable.(unique(x.ids)) +function MP.variables(p::Polynomial{T, Term{T, Monomial}}) where {T} + v = Variable[] + for t in MP.terms(p) + append!(v, MP.variables(MP.monomial(t))) + end + sort!(v, rev = true) + unique!(v) + return v +end +function MP.exponents(x::Monomial) + exps = Int[] + cur_id = nothing + for i in x.ids + if !isnothing(cur_id) && cur_id == i + exps[end] += 1 + else + cur_id = i + push!(exps, 1) + end + end + return exps +end -function Base.gcd(x::Monomial, y::Monomial) - g = Monomial(similar(x.ids, 0)) - i = j = 1 - n0 = length(x.ids) - n1 = length(y.ids) - n = n0 + n1 + 2 - T = eltype(x.ids) - M = typemax(T) - while i + j < n - xk = i <= n0 ? x.ids[i] : M - yk = j <= n1 ? y.ids[j] : M - if xk < yk - i += 1 - elseif xk > yk - j += 1 +function _degree(m::Monomial, offset) + k = 1 + while (offset + k) in eachindex(m.ids) && m.ids[offset] == m.ids[offset + k] + k += 1 + end + return k +end + +function MP.map_exponents(op, m1::Monomial, m2::Monomial) + g = Monomial(similar(m1.ids, 0)) + i = firstindex(m1.ids) + j = firstindex(m2.ids) + while i in eachindex(m1.ids) || j in eachindex(m2.ids) + if i in eachindex(m1.ids) && j in eachindex(m2.ids) && m1.ids[i] == m2.ids[j] + v = m1.ids[i] + di = _degree(m1, i) + i += di + dj = _degree(m2, j) + j += dj + elseif i in eachindex(m1.ids) && (!(j in eachindex(m2.ids)) || (m1.ids[i] < m2.ids[j])) + v = m1.ids[i] + di = _degree(m1, i) + dj = 0 + i += di else - push!(g.ids, xk) - i += 1 - j += 1 + v = m2.ids[j] + di = 0 + dj = _degree(m2, j) + j += dj + end + for _ in 1:op(di, dj) + push!(g.ids, v) end end return g end + +function MP.map_exponents!(op, m1::Monomial, m2::Monomial) + return MP.map_exponents(op, m1, m2) +end + +function MP.divides(m1::Monomial, m2::Monomial) + i = firstindex(m1.ids) + j = firstindex(m2.ids) + while i in eachindex(m1.ids) || j in eachindex(m2.ids) + if i in eachindex(m1.ids) && j in eachindex(m2.ids) && m1.ids[i] == m2.ids[j] + v = m1.ids[i] + di = _degree(m1, i) + i += di + dj = _degree(m2, j) + j += dj + elseif i in eachindex(m1.ids) && (!(j in eachindex(m2.ids)) || (m1.ids[i] < m2.ids[j])) + v = m1.ids[i] + di = _degree(m1, i) + dj = 0 + i += di + else + v = m2.ids[j] + di = 0 + dj = _degree(m2, j) + j += dj + end + if di > dj + return false + end + end + return true +end + + +function MA.promote_operation( + ::typeof(MP.substitute), + ::Type{MP.Subs}, + ::Type{Monomial}, + ::Type{Pair{Variable,T}}, +) where {T} + U = MA.promote_operation(^, T, Int) + return Term{U,Monomial} +end diff --git a/src/mpoly.jl b/src/mpoly.jl deleted file mode 100644 index c0b9fe1..0000000 --- a/src/mpoly.jl +++ /dev/null @@ -1,291 +0,0 @@ -# default polynomial type -#= -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)) -MPoly{T}(x::AbstractTerm) where T = MPoly([x]) -MPoly{T}(x::M) where {T,M<:AbstractMonomial} = MPoly(Term(x)) -MPoly{T}(x::CoeffType) where {T} = MPoly(eltype(T)(x)) -terms(x::MPoly) = x.terms -function nvariables(p::MPoly) - if monomialtype(p) <: PackedMonomial - nvariables(monomial(lt(p))) - else - error() - end -end -coeffs(x::MPoly) = (coeff(t) for t in terms(x)) - -Base.copy(x::MPoly) = MPoly(map(copy, terms(x))) -Base.one(::Type{<:MPoly{T}}) where T = MPoly(eltype(T)(1)) - -Base.:*(x::AbstractTerm, p::MPoly) = p * x -function Base.:*(p::MPoly, x::T) where {T<:AbstractTerm} - if iszero(p) - return p - elseif iszero(x) - return MPoly(emptyterm(T)) - elseif isone(x) - return p - else - MPoly(T[t * x for t in terms(p) if !iszero(t)]) - end -end -function Base.:*(p::MPoly, x::MPoly) - ts = terms(x) - out = similar(terms(p), length(ts) * length(terms(p))) - current_len = 0 - s = MPoly(out) - for tp in terms(p) - start = 1 - for tx in terms(x) - _, (current_len, start) = _add!(s, tp * tx, current_len, start) - end - end - resize!(out, current_len) - (debugmode() && !issorted(terms(s), rev=true)) && throw("Polynomial not sorted!") - - return s -end -Base.:+(p::MPoly, x::AbstractTerm) = add!(copy(p), x) -Base.:-(p::MPoly, x::AbstractTerm) = sub!(copy(p), x) - -Base.:-(p::MPoly) = -1 * p -sub!(p::MPoly, x::AbstractTerm) = add!(p, -x) - -addcoef(x::T, c) where {T<:AbstractTerm} = (c += coeff(x); return iszero(c), T(c, monomial(x))) -addcoef(x::T, c::T) where {T<:AbstractTerm} = addcoef(x, c.coeff) -subcoef(x::T, c) where {T<:AbstractTerm} = (c = coeff(x) - c; return iszero(c), T(c, monomial(x))) -subcoef(x::T, c::T) where {T<:AbstractTerm} = subcoef(x, c.coeff) - -function add!(p::MPoly, x::AbstractTerm) - q, (current_len, _) = _add!(p, x, length(terms(p)), 1); - resize!(terms(q), current_len) - (debugmode() && !issorted(terms(q), rev=true)) && throw("Polynomial not sorted!") - q -end - -@inline function _add!(p::MPoly, x::AbstractTerm, current_len, start) - iszero(x) && return p, (current_len, 1) - ts = terms(p) - i = start - for i in start:current_len - t = ts[i] - if t < x - current_len += 1 - if current_len > length(ts) - insert!(ts, i, x) - else - @inbounds @simd for j in current_len:-1:i+1 - ts[j] = ts[j-1] - end - ts[i] = x - end - return p, (current_len, i) - elseif ismatch(t, x) - return p, _add_term_matched(ts, i, t, x, current_len) - end - end - current_len += 1 - if current_len > length(ts) - push!(ts, x) - else - ts[current_len] = x - end - return p, (current_len, i) -end - -@inline function _add_term_matched(ts, i, t, x, current_len) - iz, t = addcoef(t, x) - if iz - # deleteat is somehow faster - deleteat!(ts, i) - return current_len-1, max(1, i-1) - #current_len -= 1 - #@inbounds @simd for j in i:current_len - # ts[j] = ts[j+1] - #end - #return current_len, max(1, i-1) - else - ts[i] = t - return current_len, i - end -end - -@inline function _add_rev!(p::MPoly, x::AbstractTerm, _end=1) - iszero(x) && return p, 1 - ts = terms(p) - i = _end - for i in length(ts):-1:_end - t = ts[i] - if t > x - insert!(ts, i+1, x) - return p, i - elseif ismatch(t, x) - return p, _add_term_matched_rev(ts, i, t, x) - end - end - push!(ts, x) - return p, i -end - -function _add_term_matched_rev(ts, i, t, x) - iz, t = addcoef(t, x) - if iz - deleteat!(ts, i) - return max(1, i-1) - else - ts[i] = t - return i - end -end - -function add!(p::AbstractPolynomial, x::AbstractPolynomial) - tp = terms(p) - tx = terms(x) - if tp isa Vector - sizehint!(tp, length(tp) + length(tx)) - end - for t in tx - add!(p, t) - end - return p -end -function sub!(p::AbstractPolynomial, x::AbstractPolynomial) - tp = terms(p) - tx = terms(x) - if tp isa Vector - sizehint!(tp, length(tp) + length(tx)) - end - for t in tx - sub!(p, t) - end - return p -end -Base.:+(p::AbstractPolynomial, x::AbstractPolynomial) = add!(copy(p), x) -Base.:-(p::AbstractPolynomial, x::AbstractPolynomial) = sub!(copy(p), x) - -lt(p::AbstractPolynomial) = first(terms(p)) -lc(p::AbstractPolynomial) = coeff(lt(p)) -function rmlt!(p::MPoly) - ts = terms(p) - popfirst!(ts) - return p -end -function takelt!(p::MPoly, x::MPoly) - add!(p, lt(x)) - rmlt!(x) - return p -end - -function fnmadd!(p::MPoly, l, s) - for li in terms(l) - sub!(p, li * s) - end -end - -function Base.divrem(p::MPoly, d::MPoly) - p = copy(p) - q = MPoly(similar(terms(p), 0)) - r = MPoly(similar(terms(p), 0)) - while !isempty(terms(p)) - nx, fail = lt(p) / lt(d) - if fail - takelt!(r, p) - else - #p -= d * nx - #q += nx - fnmadd!(p, d, nx) - add!(q, nx) - end - end - return q, r -end - -function divexact(p::MPoly, d::MPoly) - p = copy(p) - ts = similar(terms(p), 0) - sizehint!(ts, length(terms(d))-1) - q = MPoly(ts) - _end = 1 - while !isempty(terms(p)) - nx, fail = lt(p) / lt(d) - @assert !fail - fnmadd!(p, d, nx) - _, _end = _add_rev!(q, nx, _end) - end - (debugmode() && !issorted(terms(q), rev=true)) && throw("Polynomial not sorted!") - return q -end - -function Base.rem(p::AbstractPolynomial, d::AbstractPolynomial) - p = copy(p) - r = MPoly(similar(terms(p), 0)) - while !isempty(terms(p)) - nx, fail = lt(p) / lt(d) - if fail - takelt!(r, p) - else - sub!(p, d * nx) - end - end - return r -end - -# TODO -Base.:(/)(x::MPoly, y::MPoly) = divexact(x, y) - -Base.gcd(x::Union{AbstractTerm,AbstractMonomial}, y::MPoly) = gcd(MPoly(x), y) -Base.gcd(x::MPoly, y::Union{AbstractTerm,AbstractMonomial}) = gcd(y, x) -function Base.gcd(x::MPoly, y::MPoly) - # trival case - if iszero(x) || isone(y) - return y - elseif iszero(y) || isone(x) || x == y - return x - end - - v1, p1 = to_univariate(x) - v2, p2 = to_univariate(y) - if v1 < v2 - x, y = y, x - v1, v2 = v2, v1 - p1, p2 = p2, p1 - end - # v2 < v1 - # both are constants - v2 == NOT_A_VAR && return MPoly(gcd(lt(x), lt(y))) - if v2 < v1 - # `v2` in p2 doesn't exist in `x`, so the gcd at this level is 1 and we - # just move on to the next level - return gcd(x, content(p2)) - end - v1 == v2 || error("unreachable") - - g = gcd(p1, p2) - return univariate_to_multivariate(g) -end - -function pick_var(x::MPoly) - ts = terms(x) - v = NOT_A_VAR - for (i, t) in enumerate(ts) - if degree(t) > 0 - m = monomial(t) - vv = firstid(m) # get the minvar - if vv < v - v = vv - end - end - end - return IDType(v) -end - -function to_univariate(x::MPoly) - v = pick_var(x) - v, (v == NOT_A_VAR ? nothing : SparsePoly(x, v)) -end diff --git a/src/packedmonomial.jl b/src/packedmonomial.jl index b4d4062..d535ea9 100644 --- a/src/packedmonomial.jl +++ b/src/packedmonomial.jl @@ -36,7 +36,7 @@ new_E(::Val{E}) where {E} = max(Base._nextpow2(1+E)-1, 7) Bit packed monomial with maximum of L variables and E bits of exponents. """ -struct PackedMonomial{L,E,K} <: AbstractMonomial +struct PackedMonomial{L,E,K} <: MP.AbstractMonomial bits::NTuple{K,UInt64} function PackedMonomial{L,E}() where {L,E} EN = new_E(Val(E)) @@ -50,8 +50,8 @@ struct PackedMonomial{L,E,K} <: AbstractMonomial 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 MultivariatePolynomials.exponents(m::PackedMonomial{L,E,K}) where {L,E,K} +nvariables(::PackedMonomial{L}) where L = L +function MP.exponents(m::PackedMonomial{L,E,K}) where {L,E,K} tup = m.bits k = 0 vpu = var_per_UInt64(Val(E)) @@ -149,7 +149,7 @@ 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} +function MP.compare(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 @@ -184,7 +184,7 @@ function Base.:*(x::T, y::T) where {L,E,T<:PackedMonomial{L,E}} return T(xys) end -function MultivariatePolynomials.divides(y::T, x::T) where {L,E,T<:PackedMonomial{L,E}} +function MP.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)) @@ -193,23 +193,23 @@ function MultivariatePolynomials.divides(y::T, x::T) where {L,E,T<:PackedMonomia return o == zero(UInt64) end -MultivariatePolynomials.constantmonomial(::M) where {M<:PackedMonomial} = constantmonomial(M) -function MultivariatePolynomials.constantmonomial(M::Type{<:PackedMonomial}) +MP.constant_monomial(::M) where {M<:PackedMonomial} = MP.constant_monomial(M) +function MP.constant_monomial(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) +function MP.map_exponents!(f::F, m1::M, m2::M) where {F<:Function, L, E, M<:PackedMonomial{L,E}} + MP.map_exponents(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) +function MP.map_exponents(f::F, m1::M, m2::M) where {F<:Function, L, E, M<:PackedMonomial{L,E}} + e1 = MP.exponents(m1) + e2 = MP.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}} +function MP.div_multiple(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)) @@ -274,56 +274,36 @@ function rmid(x::T, id) where {L,E,K,T<:PackedMonomial{L,E,K}} return T(bits) end -function Base.show(io::IO, m::PackedMonomial{L,E}) where {L,E} - tup = m.bits - k = 0 - vpu = var_per_UInt64(Val(E)) - isone(m) && (print(io, '1'); return) - for i in eachindex(tup) - int = tup[i] << ((i == 1) * (E+1)) - for j in 1:vpu - (i == 1) - exponent = int >> ((E+1)*(vpu-1)) - if exponent > 0 - print(io, 'x') - print(io, int2lowerscript(k)) - exponent > 1 && print(io, int2superscript(exponent)) - end - k += 1 - int <<= (E+1) - L == k && break - end - end -end - -struct Variable{L,E} <: AbstractVariable +struct PackedVariable{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} +MP.variable_union_type(::Type{<:PackedMonomial{L,E}}) where {L,E} = PackedVariable{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) +MP.variables(::Type{<:PackedMonomial{L,E}}) where {L, E} = ntuple(i->PackedVariable{L,E}(i-1), Val(L)) +MP.variables(::M) where {M<:PackedMonomial} = MP.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}}} +MP.term_type(M::Type{<:Union{Monomial,PackedMonomial}}, ::Type{T}) where {T} = Term{T, M} +MP.polynomial_type(::Type{Term{T, M}}) where {T, M<:Union{Monomial,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 +MP.nvariables(::Polynomial{T, Term{T, M}}) where {T, L, M<:PackedMonomial{L}} = L +MP.variables(::Polynomial{T, Term{T, M}}) where {T, M<:PackedMonomial} = MP.variables(M) +Base.:(^)(x::PackedVariable{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 +function MP.monomial_type(::Type{PackedVariable{L,E}}) where {L,E} + EN = new_E(Val(E)) + K = calc_K(Val(L),Val(EN)) + return PackedMonomial{L,EN,K} +end +function MA.promote_operation( + ::typeof(MP.substitute), + ::Type{MP.Subs}, + ::Type{PackedMonomial{L,E,K}}, + ::Type{Pair{PackedVariable{L,E},T}}, +) where {L,E,K,T} + U = MA.promote_operation(^, T, Int) + return Term{U,PackedMonomial{L,E,K}} end diff --git a/src/poly.jl b/src/poly.jl deleted file mode 100644 index 72c6889..0000000 --- a/src/poly.jl +++ /dev/null @@ -1,435 +0,0 @@ -struct Uninomial <: AbstractMonomial - v::IDType - d::UInt -end - -Base.isless(x::Uninomial, y::Uninomial) = isless(degree(x), degree(y)) -degree(m::Uninomial) = Int(m.d) -var(m::Uninomial) = m.v - -coeff(m::Uninomial) = 1 -Base.one(m::Uninomial) = Uninomial(m.v, 0) -Base.show(io::IO, m::Uninomial) = print_single_monomial(io, var(m), degree(m)) - -struct Uniterm{T,M<:AbstractMonomial} <: AbstractTerm - coeff::T - uninomial::M -end -Base.convert(::Type{<:Uniterm}, t::Uninomial) = Uniterm(t) -#Base.convert(::Type{<:Uniterm}, t::CoeffType) = Uniterm(t, Uninomial()) -Uniterm(t::Uninomial) = Uniterm(coeff(t), t) - -coeff(t::Uniterm) = t.coeff -monomial(t::Uniterm) = t.uninomial -var(t::Uniterm) = var(monomial(t)) - -struct SparsePoly{T} <: AbstractPolynomial - coeffs::Vector{T} - exps::Vector{UInt} - v::IDType - - SparsePoly(coeffs::Vector{T}, exps::Vector{UInt}, v::IDType) where T = SparsePoly{T}(coeffs, exps, v) - function SparsePoly{T}(coeffs::Vector{T}, exps::Vector{UInt}, v::IDType) where T - length(coeffs) == length(exps) || error("coeffs and exps' length must match!") - new{T}(coeffs, exps, v) - end -end -Base.convert(::Type{<:SparsePoly}, m::Uninomial) = SparsePoly(Uniterm(m)) -Base.convert(::Type{<:SparsePoly}, t::Uniterm) = SparsePoly(t) -#Base.convert(::Type{<:SparsePoly}, t::CoeffType) = SparsePoly(convert(Uniterm, t)) -SparsePoly(t::Uniterm) = SparsePoly([coeff(t)], UInt[degree(t)], var(t)) -SparsePoly(c::Number, v) = SparsePoly([c], [zero(UInt)], v) -const EMPTY_EXPS = UInt[] -#SparsePoly(t::Uniterm, id::IDType) = SparsePoly(typeof(coeff(t))[], EMPTY_EXPS, id) -Base.similar(p::SparsePoly) = SparsePoly(similar(coeffs(p)), similar(p.exps), var(p)) - -var(p::SparsePoly) = p.v -coeffs(p::SparsePoly) = p.coeffs -coeff(p::SparsePoly, i) = coeffs(p)[i] -term(p::SparsePoly, i) = Uniterm(p.coeffs[i], Uninomial(var(p), p.exps[i])) -terms(p::SparsePoly) = (term(p, i) for i in eachindex(coeffs(p))) -lt(p::SparsePoly) = term(p, 1) -lc(p::SparsePoly) = coeff(p, 1) -degree(p::SparsePoly) = iszero(p) ? -1 : Int(p.exps[1]) -Base.iszero(p::SparsePoly) = isempty(p.exps) || iszero(lt(p)) - -Base.zero(t::SparsePoly) = SparsePoly(empty(t.coeffs), empty(t.exps), var(t)) -# TODO: this fails when t is zero -Base.one(t::SparsePoly) = one(lt(t)) - -Base.deleteat!(p::SparsePoly, i::Int) = (deleteat!(p.coeffs, i); deleteat!(p.exps, i); nothing) -Base.copy(p::SparsePoly) = SparsePoly(copy(coeffs(p)), copy(p.exps), var(p)) -Base.copy(p::SparsePoly{<:AbstractPolynomial}) = SparsePoly(map(copy, coeffs(p)), copy(p.exps), var(p)) -Base.:(==)(p::SparsePoly, q::SparsePoly) = all(x->x[1] == x[2], zip(terms(p), terms(q))) - -function check_poly(x, y) - v = var(x) - v == var(y) || error("$x and $y contain different variables!") - nothing -end - -function print_poly_term(io, t) - need_par = !isone(monomial(t)) - need_par && print(io, '(') - show(io, coeff(t)) - need_par && print(io, ')') - need_par && show(io, monomial(t)) - nothing -end - -function Base.show(io::IO, p::SparsePoly{<:AbstractPolynomial}) - ts = terms(p) - if isempty(ts) - print(io, 0) - return - end - n = length(ts) - t1, tr = Iterators.peel(ts) - print_poly_term(io, t1) - for t in tr - print(io, " + ") - print_poly_term(io, t) - end -end - -######################### -# Arithmetic/Algorithms # -######################### - -Base.:(^)(x::Uninomial, n::Integer) = Uninomial(var(x), degree(x)+n-1) - -function Base.:(+)(x::Uninomial, y::Uninomial) - check_poly(x, y) - if x < y - x, y = y, x - end - SparsePoly([1, 1], [degree(x), degree(y)], var(x)) -end -function Base.:(+)(x::Uniterm, y::Uniterm) - check_poly(x, y) - if ismatch(x, y) - c = x.coeff + y.coeff - return iszero(c) ? SparsePoly(typeof(x)[], EMPTY_EXPS, var(x)) : SparsePoly(Uniterm(c, monomial(x))) - else - if iszero(x) && iszero(y) - c = coeff(y) - return SparsePoly(typeof(c)[], EMPTY_EXPS, var(x)) - elseif iszero(x) - exp = degree(y) - c = coeff(y) - return SparsePoly([c], UInt[exp], var(x)) - elseif iszero(y) - exp = degree(x) - c = coeff(x) - return SparsePoly([c], UInt[exp], var(x)) - else - if x < y - x, y = y, x - end - return SparsePoly([coeff(x), coeff(y)], UInt[degree(x), degree(y)], var(x)) - end - end -end -Base.:(*)(x::CoeffType, y::Uninomial) = Uniterm(x, y) -Base.:(*)(x::Uninomial, y::CoeffType) = y * x - -function Base.:(*)(x::Uninomial, y::Uninomial) - check_poly(x, y) - Uninomial(var(x), degree(x) + degree(y)) -end -Base.:(*)(x::Uniterm, y::Uniterm) = Uniterm(coeff(x) * coeff(y), monomial(x) * monomial(y)) - -function Base.:*(p::SparsePoly, x::Uniterm) - check_poly(p, x) - if iszero(x) - return SparsePoly(x) - elseif isone(x) - return p - else - cfs = similar(coeffs(p)) - exps = similar(p.exps) - for (i, t) in enumerate(terms(p)) - nt = t * x - cfs[i] = coeff(nt) - exps[i] = degree(nt) - end - return SparsePoly(cfs, exps, var(x)) - end -end -Base.:*(x::CoeffType, y::Uniterm) = Uniterm(x * coeff(y), monomial(y)) -Base.:*(x::Uniterm, y::CoeffType) = y * x - -function Base.:*(x::SparsePoly, y::SparsePoly) - if iszero(x) - return x - elseif iszero(y) - return y - else - return sum(t->x * t, terms(y)) - end -end -smul(x, y::SparsePoly) = SparsePoly(x * coeffs(y), y.exps, var(y)) -Base.:*(x::CoeffType, y::SparsePoly) = smul(x, y) -Base.:*(x::SparsePoly, y::CoeffType) = y * x -Base.:*(x::T, y::SparsePoly{T}) where {T<:AbstractPolynomial} = smul(x, y) -Base.:*(x::SparsePoly{T}, y::T) where {T<:AbstractPolynomial} = y * x - -addcoef(x::Uniterm, c) = (c += coeff(x); return iszero(c), Uniterm(c, monomial(x))) -addcoef(x::Uniterm, c::Uniterm) = addcoef(x, coeff(c)) -Base.:(-)(x::Uniterm) = Uniterm(-coeff(x), monomial(x)) -sub!(p::SparsePoly, x::Uniterm) = add!(p, -x) -function add!(p::SparsePoly, x::Uniterm) - iszero(x) && return p - for (i, t) in enumerate(terms(p)) - if ismatch(t, x) - iz, t = addcoef(t, x) - if iz - deleteat!(p, i) - else - p.coeffs[i] = coeff(t) - end - return p - elseif t < x - insert!(p.coeffs, i, coeff(x)) - insert!(p.exps, i, degree(x)) - return p - end - end - push!(p.coeffs, coeff(x)) - push!(p.exps, degree(x)) - return p -end - -Base.div(x::Uninomial, y::Uninomial) = (x/y)[1] -function Base.:(/)(x::Uninomial, y::Uninomial) - check_poly(x, y) - xd, yd = degree(x), degree(y) - xd >= yd ? (Uninomial(var(x), xd - yd), false) : (x, true) -end - -# dest = (p * a).^n -function mulpow!(dest, p::SparsePoly, a, n::Integer) - @assert n >= 0 - for i in eachindex(p.exps) - dest.coeffs[i] = p.coeffs[i] * a - dest.exps[i] = p.exps[i] + n - end - return dest -end - -function mul!(p::SparsePoly, l) - cs = coeffs(p) - for i in eachindex(cs) - cs[i] *= l - end -end - -function pseudorem(p::SparsePoly, d::SparsePoly) - check_poly(p, d) - degree(p) < degree(d) && return p - k = degree(p) - degree(d) + 1 - l = lc(d) - dd = similar(d) - p = copy(p) - while !iszero(p) && degree(p) >= degree(d) - s = mulpow!(dd, d, lc(p), degree(p) - degree(d)) - mul!(p, l) - sub!(p, s) - k -= 1 - end - return l^k * p -end - -function termwise_content(a::MPoly) - ts = terms(a) - length(ts) == 1 && return ts[1] - g = gcd(ts[1], ts[2]) - isone(g) || for i in 3:length(ts) - g = gcd(g, ts[i]) - isone(g) && break - end - g -end - -function content(a::SparsePoly) - cfs = coeffs(a) - length(cfs) == 1 && return first(cfs) - # If any coeff here is a multivariate polynomial with only one term, then we - # just need to compute the termwize content. - MP = eltype(cfs) - if MP <: AbstractPolynomial - for i in eachindex(cfs) - ts = terms(cfs[i]) - if length(ts) == 1 - g = gcd(termwise_content(cfs[1]), termwise_content(cfs[2])) - isone(g) || for i in 3:length(cfs) - g = gcd(g, termwise_content(cfs[i])) - isone(g) && break - end - return MP(g) - end - end - end - - g = gcd(cfs[1], cfs[2]) - isone(g) || for i in 3:length(cfs) - g = gcd(g, cfs[i]) - isone(g) && break - end - g -end - -#function univariate_gcd(x::AbstractPolynomial, y::AbstractPolynomial) -# while !iszero(y) -# x = pseudorem(x, y) -# x, y = y, x -# end -# return x#, a / x, b / x -#end - -function Base.divrem(p::SparsePoly{T}, d::SparsePoly{T}) where T<:CoeffType - p = copy(p) - q = zero(p) - r = zero(p) - while !iszero(p) - nx, fail = lt(p) / lt(d) - if fail - break - else - p -= d * nx - q += nx - end - end - return q, r -end - -function divexact(a::SparsePoly{T}, b::T) where {T<:AbstractPolynomial} - cfs = [divexact(c, b) for c in coeffs(a)] - SparsePoly(cfs, a.exps, var(a)) -end - -divexact(c, b) = ((d, r) = divrem(c, b); @assert iszero(r); d;) -function divexact(a::SparsePoly{T}, b::T) where T - cfs = [divexact(c, b) for c in coeffs(a)] - SparsePoly(cfs, a.exps, var(a)) -end - -function Base.gcd(x::SparsePoly, y::SparsePoly) - if degree(y) > degree(x) - x, y = y, x - end - if iszero(y) - return x - elseif isone(y) - return y - end - c1, x = contprim(x) - c2, y = contprim(y) - c = gcd(c1, c2) - g = one(eltype(coeffs(x))) - h = one(eltype(coeffs(x))) - while true - r = pseudorem(x, y) - iszero(r) && break - degree(r) == 0 && return SparsePoly(c, var(x)) - - d = degree(x) - degree(y) - x, y = y, divexact(r, g*h^d) - g = lc(x) - h = d > 1 ? divexact(g^d, h^(d - 1)) : h^(1 - d)*g^d - end - return c*primpart(y) -end - -primpart(p) = divexact(p, content(p)) -contprim(p) = (c = content(p); (c, divexact(p, c))) - -############## -# Conversion # -############## - -# mpoly2poly(5x⁴ + 2x³ + x² + xy + y^2, 0x00001 [y]) -# -> -# ---------------- constants -# | -# y^2 + xy + (5x⁴ + 2x³ + x²) -# -# M2P(a, v) -function emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree) - v = var(poly) - if olddegree > 0 - t = ts[perm[chunk_start_idx]] - coeff = MPoly(term2polycoeff(t, v)) - for i in chunk_start_idx+1:idx-1 - t = ts[perm[i]] - add!(coeff, term2polycoeff(t, v)) - end - else - coeff = MPoly(ts[perm[chunk_start_idx]]) - for i in chunk_start_idx+1:idx-1 - add!(coeff, ts[perm[i]]) - end - end - push!(poly.exps, olddegree) - push!(poly.coeffs, coeff) - nothing -end -term2polycoeff(t::Term, v::IDType) = Term(coeff(t), rmid(monomial(t), v)) -function SparsePoly(p::MPoly, v::IDType) - ts = terms(p) - pows = map(ts) do t - degree(monomial(t), v) - end - perm = sortperm(pows, rev=true) - coeffs = typeof(p)[] - exps = UInt[] - poly = SparsePoly(coeffs, exps, v) - - olddegree = pows[perm[1]] - chunk_start_idx = idx = 1 - while idx <= length(perm) - degree = pows[perm[idx]] - if olddegree != degree # new chunk - emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree) - chunk_start_idx = idx - olddegree = degree - end - idx += 1 - end - if !isempty(chunk_start_idx:idx-1) # remember to handle the last one - emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree) - end - return poly -end - -univariate_to_multivariate(g::SparsePoly{<:AbstractPolynomial}) = univariate_to_multivariate(g, monomialtype(lc(g))) -function univariate_to_multivariate(g, ::Type{<:Monomial}) - cfs = coeffs(g) - eps = g.exps - v = var(g) - @assert !isempty(eps) - #sum(zip(cfs, eps)) do (c, e) - # c * Monomial(fill(v, e)) - #end - s = cfs[1] * Monomial(fill(v, eps[1])) - for i in 2:length(cfs) - c = cfs[i] - e = eps[i] - add!(s, c * Monomial(fill(v, e))) - end - return s -end - -function univariate_to_multivariate(g, P::Type{<:PackedMonomial}) - cfs = coeffs(g) - eps = g.exps - v = var(g) - @assert !isempty(eps) - s = cfs[1] * P(v)^eps[1] - for i in 2:length(cfs) - c = cfs[i] - e = eps[i] - add!(s, c * P(v)^e) - end - return s -end diff --git a/src/variable.jl b/src/variable.jl new file mode 100644 index 0000000..8afb2cb --- /dev/null +++ b/src/variable.jl @@ -0,0 +1,15 @@ +abstract type AbstractVariable <: MP.AbstractVariable end + +function Base.:(==)(v1::V, v2::V) where {V<:AbstractVariable} + v1 === v2 +end +function Base.isless(v1::V, v2::V) where {V<:AbstractVariable} + isless(v2.id, v1.id) +end +function MP.name_base_indices(v::AbstractVariable) + "x", (v.id,) +end +MP.name(v::AbstractVariable) = string("x", MP.unicode_subscript(v.id)) +function MP.substitute(::MP.Subs, v::V, p::Pair{V,Int64}) where {V<:AbstractVariable} + v == p[1] ? p[2] : v +end diff --git a/test/runtests.jl b/test/runtests.jl index d084e0d..0b2dba9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,32 +1,11 @@ using SIMDPolynomials -using SIMDPolynomials: divexact, lc, var +import MultivariatePolynomials as MP using Test SIMDPolynomials.debugmode() = true -@testset "Fuzz SparsePoly" begin - x = Uninomial(0, 1) - for i in 1:10000 - p = sum(rand(-2:2)*x^i for i in 0:10) - q = sum(rand(-2:2)*x^i for i in 0:5) - try - g = gcd(p, q) - pdg = divexact(p, g) - qdg = divexact(q, g) - @test gcd(pdg, qdg) == one(x) - pq = p * q - @test divexact(pq, p) == q - @test divexact(pq, q) == p - catch - display(p) - display(q) - rethrow() - end - end -end - @testset "Fuzz MPoly" begin - for monos in [[PackedMonomial{4,8}(i) for i in 0:2], [Monomial([i]) for i in 0:2]] + for monos in ([PackedMonomial{4,8}(i) for i in 0:2], [Monomial([i]) for i in 0:2]) x, y, z = monos for _ in 1:10000 pterms = [rand(-2:2)*prod(rand([x, y, z])^i for i in 0:3) for j in 0:7] @@ -34,34 +13,26 @@ end q = sum(rand(-2:2)*prod(rand([x, y, z])^i for i in 0:2) for j in 0:5) try g = gcd(p, q) - pdg = divexact(p, g) - qdg = divexact(q, g) - @test gcd(pdg, qdg) == one(x) + pdg = MP.div_multiple(p, g) + qdg = MP.div_multiple(q, g) + if iszero(p) && iszero(q) + @test gcd(pdg, qdg) == zero(x) + else + @test gcd(pdg, qdg) == one(x) + end catch - display(p) - display(q) + println(p) + println(q) rethrow() end end end end -@testset "pseudorem" begin - x = Uninomial(0, 1) - p = 2x^10 + x^7 + 7*x^2 + x + 3x - q = (p+one(p)) * (p+2one(p)) * (p+3one(p)) - @test pseudorem(q, p) == 12582912*one(p) - q = x^7 + 20*one(x) - @test pseudorem(q, p) == q - @test pseudorem(p, q) == -40*x^3 + 7*x^2 + 4*x - 20*one(x) - q = x^6 + 23*one(x) - @test pseudorem(p, q) == -46*x^4 + 7*x^2 - 19*x -end - function test_gcd(x, y) g1 = gcd(x, y) g2 = gcd(y, x) - @test sign(lc(g1)) * g1 == sign(lc(g2)) * g2 + @test sign(MP.leading_coefficient(g1)) * g1 == sign(MP.leading_coefficient(g2)) * g2 return g1 end @@ -76,12 +47,8 @@ end e3 = 7 e4 = 10 p = c1 * y^e1 + c2 * y^e2 + c3 * y^e3 + c4 * y^e4 - pp = SIMDPolynomials.SparsePoly(p, y.ids[1]) - @test var(pp) == y.ids[1] - @test coeffs(pp) == [c4, c3, c2, c1] - @test pp.exps == [e4, e3, e2, e1] q = prod(i->p + i, 0:3) - @test length(terms(q)) == 262 + @test length(MP.terms(q)) == 262 for i in 0:3 @test test_gcd(p + i, q) == p + i end @@ -126,7 +93,7 @@ end e4 = 10 p = c1 * y^e1 + c2 * y^e2 + c3 * y^e3 + c4 * y^e4 q = prod(i->p + i, 0:3); - @test length(terms(q)) == 262 + @test length(MP.terms(q)) == 262 for i in 0:3 @test test_gcd(p + i, q) == p + i end @@ -152,7 +119,7 @@ end @test test_gcd(p, p1) == p1 k = p; for i in 1:n - k = divexact(k, p1) + k = div(k, p1) end @test k == 1 q = (p + 1) * p;