Performance

.github/workflows/CI.yml

@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1.7'
           - 'nightly'
         os:

.gitignore

@@ -3,3 +3,4 @@
 *.jl.mem
 /docs/build/
 statprof/
+profile.pb.gz

Project.toml

@@ -1,12 +1,13 @@
 name = "LoopPoly"
 uuid = "4edc584b-a88b-4acf-80bb-891198a11e01"
-authors = ["Yingbo Ma <mayingbo5@gmail.com> and contributors"]
+authors = ["Yingbo Ma <mayingbo5@gmail.com> and Chris Elrod <elrodc@gmail.com>"]
 version = "0.1.0"
 
 [deps]
 SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 
 [compat]
+SIMD = "3"
 julia = "1"
 
 [extras]

src/LoopPoly.jl

@@ -7,8 +7,21 @@ export mpoly2poly
 export univariate_gcd, pseudorem
 export PackedMonomial
 
-include("mpoly.jl")
+struct Ret{V} <: Function
+    value::V
+    Ret{V}(value) where {V} = new{V}(value)
+    Ret(value) = new{Core.Typeof(value)}(value)
+end
+
+(obj::Ret)(args...; kw...) = obj.value
+
+debugmode() = false
+
+include("interface.jl")
+include("utils.jl")
+include("monomial.jl")
 include("packedmonomial.jl")
+include("mpoly.jl")
 include("poly.jl")
 
 end

src/interface.jl

@@ -0,0 +1,160 @@
+###
+### 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)
+
+###
+### 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))
+
+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[]
+
+###
+### AbstractPolynomial: terms, copy
+###
+abstract type AbstractPolynomial <: Number end
+
+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))
+
+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

src/monomial.jl

@@ -0,0 +1,148 @@
+const IDType = UInt32
+const NOT_A_VAR = typemax(IDType)
+const EMPTY_IDS = IDType[]
+
+struct Monomial <: AbstractMonomial
+    ids::Vector{IDType}
+end
+Monomial() = Monomial(EMPTY_IDS)
+degree(x::Monomial) = length(x.ids)
+Base.copy(x::Monomial) = Monomial(copy(x.ids))
+
+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)
+    ids = y.ids
+    i = j = 1
+    n0 = length(ids)
+    n1 = length(x.ids)
+    r = Monomial(similar(ids, n0 + n1))
+    T = eltype(ids)
+    M = typemax(T)
+    for k in 1:(n0 + n1)
+        a = i <= n0 ? ids[i] : M
+        b = j <= n1 ? x.ids[j] : M
+        if a < b
+            i += 1
+            r.ids[k] = a
+        else
+            j += 1
+            r.ids[k] = b
+        end
+    end
+    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
+function Base.isless(x::Monomial, y::Monomial)
+    dx = degree(x)
+    dy = degree(y)
+    dx < dy && return true
+    dx > dy && return false
+    for i in 1:dx
+        vx = x.ids[i]
+        vy = y.ids[i]
+        vx != vy && return vx > vy
+    end
+    return false
+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
+        else
+            push!(g.ids, xk)
+            i += 1
+            j += 1
+        end
+    end
+    return g
+end