diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index e062b22..510c7f1 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1.7'
           - 'nightly'
         os:
diff --git a/.gitignore b/.gitignore
index 1bf7c5d..46bfd1a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@
 *.jl.mem
 /docs/build/
 statprof/
+profile.pb.gz
diff --git a/Project.toml b/Project.toml
index 101009d..247f1b1 100644
--- a/Project.toml
+++ b/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]
diff --git a/src/LoopPoly.jl b/src/LoopPoly.jl
index 92afd81..921cc26 100644
--- a/src/LoopPoly.jl
+++ b/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
diff --git a/src/interface.jl b/src/interface.jl
new file mode 100644
index 0000000..b834bbd
--- /dev/null
+++ b/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
diff --git a/src/monomial.jl b/src/monomial.jl
new file mode 100644
index 0000000..84335c6
--- /dev/null
+++ b/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
diff --git a/src/mpoly.jl b/src/mpoly.jl
index 3525696..3feb88a 100644
--- a/src/mpoly.jl
+++ b/src/mpoly.jl
@@ -1,282 +1,58 @@
-debugmode() = false
-const Rat = Union{Rational{<:Integer},Integer}
-
-const IDType = UInt32
-const NOT_A_VAR = typemax(IDType)
-const EMPTY_IDS = IDType[]
-
-abstract type AbstractMonomial <: Number end
-
-struct Monomial <: AbstractMonomial
-    ids::Vector{IDType}
-end
-Monomial() = Monomial(EMPTY_IDS)
-Base.promote_rule(::Type{Monomial}, ::Type{<:Rat}) = Term
-
-const VARNAME_DICT = Dict{IDType, String}(0 => "x", 1 => "y", 2 => "z")
-const PRETTY_PRINT = Ref(true)
-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
-
-Base.:-(y::Monomial) = Term(-1, y)
-Base.:+(y::Monomial, x::Monomial) = Term(y) + Term(x)
-Base.:^(y::Monomial, n::Integer) = iszero(n) ? Monomial() : Base.power_by_squaring(y, n)
-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)
-
-Base.isone(x::Monomial) = iszero(degree(x))
-degree(x::Monomial) = length(x.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
-
-Base.:*(c::Rat, m::Monomial) = Term(c, m)
-Base.:*(m::Monomial, c::Rat) = c * m
-
-abstract type AbstractTerm <: Number end
-struct Term <: AbstractTerm
-    coeff::Rational{Int}
-    monomial::Monomial
-end
-monomial(x::Term) = x.monomial
-Term(x) = Term(x, Monomial())
-Term(m::Monomial) = Term(1, m)
-coeff(x::Term) = x.coeff
-degree(x::Term) = degree(monomial(x))
-
-Base.promote_rule(t::Type{<:AbstractTerm}, ::Type{<:AbstractMonomial}) = t
-Base.promote_rule(t::Type{<:AbstractTerm}, ::Type{<:Rat}) = t
-Base.:(==)(x::AbstractTerm, y::AbstractTerm) = x === y || (coeff(x) == coeff(y) && monomial(x) == monomial(y))
-
-print_coeff(io::IO, coeff) = isinteger(coeff) ? print(io, Integer(coeff)) : print(io, coeff)
-function Base.show(io::IO, x::AbstractTerm)
-    printed = false
-    if x.coeff != 1
-        if x.coeff == -1
-            print(io, '-')
-        else
-            print_coeff(io, x.coeff)
-            printed = true
-            PRETTY_PRINT[] || print(io, '*')
-        end
-    end
-    m = monomial(x)
-    if !(printed && isone(m))
-        show(io, m)
-    end
-end
-
-ismatch(x::AbstractTerm, y::AbstractTerm) = monomial(x) == monomial(y)
-Base.isless(x::Term, y::Term) = isless(monomial(x), monomial(y))
-Base.iszero(x::Term) = iszero(x.coeff)
-Base.isone(x::Term) = isone(x.coeff) && isone(monomial(x))
-addcoef(x::Term, c) = (c += x.coeff; return iszero(c), Term(c, monomial(x)))
-addcoef(x::Term, c::Term) = addcoef(x, c.coeff)
-subcoef(x::Term, c) = (c = x.coeff - c; return iszero(c), Term(c, monomial(x)))
-subcoef(x::Term, c::Term) = subcoef(x, c.coeff)
-
-Base.isinteger(x::Term) = isinteger(x.coeff)
-Base.:*(x::Term, y::Term) = Term(x.coeff * y.coeff, monomial(x) * monomial(y))
-function Base.:+(x::Term, y::Term)
-    if ismatch(x, y)
-        c = x.coeff + y.coeff
-        return iszero(c) ? MPoly(Term(0)) : MPoly(Term(c, monomial(x)))
-    else
-        return x < y ? MPoly(Term[y, x]) : MPoly(Term[x, y])
-    end
-end
-
-Base.:-(x::Term) = Term(-x.coeff, monomial(x))
-function Base.:-(x::Term, y::Term)
-    if ismatch(x, y)
-        c = x.coeff - y.coeff
-        return iszero(c) ? MPoly(Term(0)) : MPoly(Term(c, x.coeff))
-    else
-        y = -y
-        return x < y ? MPoly(Term[y, x]) : MPoly(Term[x, y])
-    end
-end
-
-function Base.:(/)(y::AbstractTerm, x::AbstractTerm)
-    m, fail = monomial(y) / monomial(x)
-    if fail
-        Term(y.coeff, m), fail
-    else
-        Term(y.coeff/x.coeff, m), fail
-    end
-end
-
-const EMPTY_TERMS = Term[]
-abstract type AbstractPoly <: Number end
-struct MPoly <: AbstractPoly
-    terms::Vector{Term}
-end
-MPoly() = MPoly(EMPTY_TERMS)
-MPoly(x::Term) = MPoly([x])
-MPoly(x::Union{Rat,Monomial}) = MPoly(Term(x))
-Base.promote_rule(p::Type{<:AbstractPoly}, ::Type{<:AbstractMonomial}) = p
-Base.promote_rule(p::Type{<:AbstractPoly}, ::Type{<:AbstractTerm}) = p
-Base.promote_rule(p::Type{<:AbstractPoly}, ::Type{<:Rat}) = p
+# default polynomial type
+struct MPoly{T} <: AbstractPolynomial
+    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
-Base.empty!(x::MPoly) = (empty!(terms(x)); x)
-function Base.show(io::IO, p::AbstractPoly)
-    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
 
-Base.copy(x::Monomial) = Monomial(copy(x.ids))
-Base.copy(x::Term) = Term(copy(coeff(x)), copy(monomial(x)))
 Base.copy(x::MPoly) = MPoly(map(copy, terms(x)))
-#function largercopy(x::MPoly, i::Int)
-#    n = length(terms(x))
-#    terms = similar(terms(x), i + n)
-#    copyto!(terms, 1, terms(x), 1, n)
-#    MPoly(terms)
-#end
-Base.:(==)(x::MPoly, y::MPoly) = x === y || (terms(x) == terms(y))
-Base.iszero(x::MPoly) = isempty(terms(x))
-Base.isone(x::MPoly) = (ts = terms(x); length(ts) == 1 && isone(only(ts)))
+Base.one(::Type{<:MPoly{T}}) where T = MPoly(eltype(T)(1))
 
-Base.:*(x::Term, p::MPoly) = p * x
-function Base.:*(p::MPoly, x::Term)
+Base.:*(x::AbstractTerm, p::MPoly) = p * x
+function Base.:*(p::MPoly, x::T) where {T<:AbstractTerm}
     if iszero(x)
-        return MPoly()
+        return MPoly(emptyterm(T))
     elseif isone(x)
         return p
     else
-        ts = Term[t * x for t in terms(p) if !iszero(t)]
-        MPoly(ts)
+        MPoly(T[t * x for t in terms(p) if !iszero(t)])
     end
 end
 function Base.:*(p::MPoly, x::MPoly)
-    sum(t->p * t, terms(x))
+    ts = terms(x)
+    out = similar(terms(p), length(ts) * length(terms(p)))
+    resize!(out, 0)
+    s = MPoly(out)
+    for tp in terms(p)
+        for tx in terms(x)
+            add!(s, tp * tx)
+        end
+    end
+
+    return s
 end
-Base.:+(p::MPoly, x::Term) = add!(copy(p), x)
-Base.:-(p::MPoly, x::Term) = sub!(copy(p), x)
+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::Term) = add!(p, -x)
+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::Term)
+function add!(p::MPoly, x::AbstractTerm)
     q = _add!(p, x);
     (debugmode() && !issorted(terms(q), rev=true)) && throw("Polynomial not sorted!")
     q
 end
-function _add!(p::MPoly, x::Term)
+function _add!(p::MPoly, x::AbstractTerm)
     iszero(x) && return p
     ts = terms(p)
     for (i, t) in enumerate(ts)
@@ -297,23 +73,23 @@ function _add!(p::MPoly, x::Term)
     return p
 end
 
-function add!(p::AbstractPoly, x::AbstractPoly)
+function add!(p::AbstractPolynomial, x::AbstractPolynomial)
     for t in terms(x)
         add!(p, t)
     end
     return p
 end
-function sub!(p::AbstractPoly, x::AbstractPoly)
+function sub!(p::AbstractPolynomial, x::AbstractPolynomial)
     for t in terms(x)
         sub!(p, t)
     end
     return p
 end
-Base.:+(p::AbstractPoly, x::AbstractPoly) = add!(copy(p), x)
-Base.:-(p::AbstractPoly, x::AbstractPoly) = sub!(copy(p), x)
+Base.:+(p::AbstractPolynomial, x::AbstractPolynomial) = add!(copy(p), x)
+Base.:-(p::AbstractPolynomial, x::AbstractPolynomial) = sub!(copy(p), x)
 
-lt(p::AbstractPoly) = first(terms(p))
-lc(p::AbstractPoly) = coeff(lt(p))
+lt(p::AbstractPolynomial) = first(terms(p))
+lc(p::AbstractPolynomial) = coeff(lt(p))
 function rmlt!(p::MPoly)
     ts = terms(p)
     popfirst!(ts)
@@ -355,7 +131,7 @@ function divexact(p::MPoly, d::MPoly)
     return q
 end
 
-function Base.rem(p::AbstractPoly, d::AbstractPoly)
+function Base.rem(p::AbstractPolynomial, d::AbstractPolynomial)
     p = copy(p)
     r = MPoly(similar(terms(p), 0))
     while !isempty(terms(p))
@@ -369,41 +145,11 @@ function Base.rem(p::AbstractPoly, d::AbstractPoly)
     return r
 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
-
-function Base.gcd(x::Term, y::Term)
-    #g, a, b = gcd(monomial(x), monomial(y))
-    g = gcd(monomial(x), monomial(y))
-    gr = gcd(x.coeff, y.coeff)
-    return Term(gr, g)#, Term(x.coeff / gr, a), Term(y.coeff / gr, b)
-end
-
+# TODO
 Base.:(/)(x::MPoly, y::MPoly) = divexact(x, y)
 
-Base.gcd(x::Term, y::MPoly) = gcd(MPoly(x), y)
-Base.gcd(x::MPoly, y::Term) = gcd(y, x)
+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)
@@ -439,13 +185,13 @@ function pick_var(x::MPoly)
     for (i, t) in enumerate(ts)
         if degree(t) > 0
             m = monomial(t)
-            vv = m.ids[1] # get the minvar
+            vv = firstid(m) # get the minvar
             if vv < v
                 v = vv
             end
         end
     end
-    return v
+    return IDType(v)
 end
 
 function to_univariate(x::MPoly)
diff --git a/src/packedmonomial.jl b/src/packedmonomial.jl
index 8a79412..379999b 100644
--- a/src/packedmonomial.jl
+++ b/src/packedmonomial.jl
@@ -34,12 +34,17 @@ struct PackedMonomial{L,E,K} <: AbstractMonomial
     bits::NTuple{K,UInt64}
     function PackedMonomial{L,E}() where {L,E}
         EN = new_E(Val(E))
-        PackedMonomial{L,EN}(ntuple(Returns(zero(UInt64)), calc_K(Val(L),Val(EN))))
+        K = calc_K(Val(L),Val(EN))
+        PackedMonomial{L,EN,K}()
+    end
+    function PackedMonomial{L,E,K}() where {L,E,K}
+        PackedMonomial{L,E,K}(ntuple(Ret(zero(UInt64)), Val(K)))
     end
     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
 
+PackedMonomial{L,E,K}(i::Integer) where {L,E,K} = PackedMonomial{L,E}(i)
 function PackedMonomial{L,OE}(i::Integer) where {L,OE} # x_i
     @assert i < L
     E = new_E(Val(OE))
@@ -58,8 +63,6 @@ function PackedMonomial{L,OE}(i::Integer) where {L,OE} # x_i
     PackedMonomial{L,E}(t)
 end
 
-
-
 function degree(p::PackedMonomial{L,E}) where {L,E}
     (p.bits[1] >> ((E+1) * (var_per_UInt64(Val(E))-1))) % Int
 end
@@ -142,63 +145,6 @@ end
     zero_bits(x, Val(E)) != zero(UInt64)
 end
 
-@generated function reduce_tup(f::F, inds::Tuple{Vararg{Any,N}}) where {F,N}
-    q = Expr(:block, Expr(:meta, :inline, :propagate_inbounds))
-    if N == 1
-        push!(q.args, :(inds[1]))
-        return q
-    end
-    syms = Vector{Symbol}(undef, N)
-    i = 0
-    for n ∈ 1:N
-        syms[n] = iₙ = Symbol(:i_, (i += 1))
-        push!(q.args, Expr(:(=), iₙ, Expr(:ref, :inds, n)))
-    end
-    W =  1 << (8sizeof(N) - 2 - leading_zeros(N))
-    while W > 0
-        _N = length(syms)
-        for _ ∈ 2W:W:_N
-            for w ∈ 1:W
-                new_sym = Symbol(:i_, (i += 1))
-                push!(q.args, Expr(:(=), new_sym, Expr(:call, :f, syms[w], syms[w+W])))
-                syms[w] = new_sym
-            end
-            deleteat!(syms, 1+W:2W)
-        end
-        W >>>= 1
-    end
-    q
-end
-
-struct GetWithShift{F,T}
-    f::F
-    tup::T
-    shift::Int
-end
-(g::GetWithShift)(i) = g.f(g.tup[i + g.shift])
-
-struct GetWithShift2{F,T}
-    f::F
-    tup1::T
-    tup2::T
-end
-(g::GetWithShift2)(i) = g.f(g.tup1[i], g.tup2[i])
-struct GetWithShift3{F,T,S}
-    f::F
-    tup::T
-    s::S
-end
-(g::GetWithShift3)(i) = g.f(g.tup[i], g.s)
-
-function _fmap(f::F, x::Tuple{Vararg{Any,K}}) where {K,F}
-    ntuple(GetWithShift(f, x, 0), Val(K))
-end
-function _fmap(f::F, x::Tuple{Vararg{Any,K}}, y::Tuple{Vararg{Any,K}}) where {K,F}
-    ntuple(GetWithShift2(f, x, y), Val(K))
-end
-function _fmap(f::F, x::Tuple{Vararg{Any,K}}, y) where {K,F}
-    ntuple(GetWithShift3(f, x, y), Val(K))
-end
 function Base.:*(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))
@@ -210,8 +156,8 @@ function Base.:/(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) && overflowed_error()
+    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)
@@ -220,6 +166,55 @@ function Base.:^(x::T, y::Integer) where {L,E,T<:PackedMonomial{L,E}}
     o != zero(UInt64) && overflowed_error()
     return T(xys)
 end
+zero_nondegreemask(::Val{7 }) = 0xff00000000000000
+zero_nondegreemask(::Val{15}) = 0xffff000000000000
+zero_nondegreemask(::Val{31}) = 0xffffffff00000000
+zero_nondegreemask(::Val{63}) = 0xffffffffffffffff
+function firstid(m::PackedMonomial{L,E,K}) where {L,E,K}
+    b = m.bits[1] & (~zero_nondegreemask(Val(E)))
+    if (b != zero(b))
+        return (leading_zeros(b) ÷ (E+1)) - 1
+    end
+    vpu = var_per_UInt64(Val(E))
+    b = vpu - 1
+    for k in 2:K
+        bk = m.bits[k]
+        if (bk != zero(bk))
+            return (leading_zeros(bk) ÷ (E+1)) + b
+        end
+        b += vpu
+    end
+    error("firstTermID should only be called if degree > 0.");
+    return 0
+end
+
+function degree(m::PackedMonomial{L,E,K}, id) where {L,E,K}
+    vpu = var_per_UInt64(Val(E))
+    d = K == 1 ? 0 : (id + 1) ÷ vpu
+    r = (id + 1) - (d * vpu)
+    b = m.bits[d+1] << (r * (E + 1))
+    return b >> ((E + 1) * (vpu - 1))
+end
+
+function rmid(x::T, id) where {L,E,K,T<:PackedMonomial{L,E,K}}
+    bits = x.bits
+    vpu = var_per_UInt64(Val(E))
+    d = K == 1 ? 0 : (id + 1) ÷ vpu
+    r = (id + 1) - (d * vpu)
+    m = zero_nondegreemask(Val(E)) >> (r*(E+1))
+    oldBits = bits[d+1]
+    b = oldBits & (~m)
+    remDegree = (oldBits & m) >> ((E+1)*(vpu-1-r))
+
+    o = remDegree << ((E + 1) * (vpu - 1))
+    if d != zero(d)
+        bits = Base.setindex(bits, b, d+1)
+        bits = Base.setindex(bits, bits[1] - o, 1)
+    else
+        bits = Base.setindex(bits, b - o, 1)
+    end
+    return T(bits)
+end
 
 function Base.show(io::IO, m::PackedMonomial{L,E}) where {L,E}
     tup = m.bits
diff --git a/src/poly.jl b/src/poly.jl
index 3663250..254a22a 100644
--- a/src/poly.jl
+++ b/src/poly.jl
@@ -1,39 +1,29 @@
-# coeff, length, shift_left
-
 struct Uninomial <: AbstractMonomial
     v::IDType
     d::UInt
 end
+
 Base.isless(x::Uninomial, y::Uninomial) = isless(degree(x), degree(y))
-degree(m::Uninomial) = m.d
+degree(m::Uninomial) = Int(m.d)
 var(m::Uninomial) = m.v
 
 coeff(m::Uninomial) = 1
-Base.isone(m::Uninomial) = iszero(degree(m))
 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} <: AbstractTerm
+struct Uniterm{T,M<:AbstractMonomial} <: AbstractTerm
     coeff::T
-    uninomial::Uninomial
+    uninomial::M
 end
 Base.convert(::Type{<:Uniterm}, t::Uninomial) = Uniterm(t)
-#Base.convert(::Type{<:Uniterm}, t::Rat) = Uniterm(t, Uninomial())
+#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))
 
-degree(t::Uniterm) = degree(monomial(t))
-Base.isless(x::Uniterm, y::Uniterm) = isless(monomial(x), monomial(y))
-
-Base.iszero(t::Uniterm) = iszero(coeff(t))
-Base.isone(x::Uniterm) = isone(coeff(x)) && isone(monomial(x))
-Base.zero(t::Uniterm) = Uniterm(zero(coeff(t)), one(monomial(t)))
-Base.one(t::Uniterm) = Uniterm(one(coeff(t)), one(monomial(t)))
-
-struct SparsePoly{T} <: AbstractPoly
+struct SparsePoly{T} <: AbstractPolynomial
     coeffs::Vector{T}
     exps::Vector{UInt}
     v::IDType
@@ -46,11 +36,11 @@ struct SparsePoly{T} <: AbstractPoly
 end
 Base.convert(::Type{<:SparsePoly}, m::Uninomial) = SparsePoly(Uniterm(m))
 Base.convert(::Type{<:SparsePoly}, t::Uniterm) = SparsePoly(t)
-#Base.convert(::Type{<:SparsePoly}, t::Rat) = SparsePoly(convert(Uniterm, t))
-SparsePoly(t::Uniterm) = SparsePoly([coeff(t)], [degree(t)], var(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)
+#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
@@ -60,7 +50,7 @@ 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 : p.exps[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) = zero(lt(t))
@@ -68,7 +58,7 @@ 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{<:AbstractPoly}) = SparsePoly(map(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)
@@ -86,7 +76,7 @@ function print_poly_term(io, t)
     nothing
 end
 
-function Base.show(io::IO, p::SparsePoly{<:AbstractPoly})
+function Base.show(io::IO, p::SparsePoly{<:AbstractPolynomial})
     ts = terms(p)
     if isempty(ts)
         print(io, 0)
@@ -118,16 +108,29 @@ function Base.:(+)(x::Uniterm, y::Uniterm)
     check_poly(x, y)
     if ismatch(x, y)
         c = x.coeff + y.coeff
-        return iszero(c) ? SparsePoly(x, var(x)) : SparsePoly(Uniterm(c, monomial(x)))
+        return iszero(c) ? SparsePoly(typeof(x)[], EMPTY_EXPS, var(x)) : SparsePoly(Uniterm(c, monomial(x)))
     else
-        if x < y
-            x, y = y, x
+        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
-        return SparsePoly([coeff(x), coeff(y)], [degree(x), degree(y)], var(x))
     end
 end
-Base.:(*)(x::Rat, y::Uninomial) = Uniterm(x, y)
-Base.:(*)(x::Uninomial, y::Rat) = y * x
+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)
@@ -152,15 +155,15 @@ function Base.:*(p::SparsePoly, x::Uniterm)
         return SparsePoly(cfs, exps, var(x))
     end
 end
-Base.:*(x::Rat, y::Uniterm) = Uniterm(x * coeff(y), monomial(y))
-Base.:*(x::Uniterm, y::Rat) = y * x
+Base.:*(x::CoeffType, y::Uniterm) = Uniterm(x * coeff(y), monomial(y))
+Base.:*(x::Uniterm, y::CoeffType) = y * x
 
 Base.:*(x::SparsePoly, y::SparsePoly) = sum(t->x * t, terms(y))
 smul(x, y::SparsePoly) = SparsePoly(x * coeffs(y), y.exps, var(y))
-Base.:*(x::Rat, y::SparsePoly) = smul(x, y)
-Base.:*(x::SparsePoly, y::Rat) = y * x
-Base.:*(x::T, y::SparsePoly{T}) where {T<:AbstractPoly} = smul(x, y)
-Base.:*(x::SparsePoly{T}, y::T) where {T<:AbstractPoly} = y * x
+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))
@@ -192,7 +195,7 @@ 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(xd - yd, var(x)), false) : (x, true)
+    xd >= yd ? (Uninomial(var(x), xd - yd), false) : (x, true)
 end
 
 # dest = (p * a).^n
@@ -205,50 +208,27 @@ function mulpow!(dest, p::SparsePoly, a, n::Integer)
     return dest
 end
 
-function fnmadd!(p::SparsePoly, l, s)
-    for li in terms(l)
-        sub!(p, li * s)
-    end
-end
+#function fnmadd!(p::SparsePoly, l, s)
+#    for li in terms(l)
+#        sub!(p, li * s)
+#    end
+#end
 function mul!(p::SparsePoly, l)
     cs = coeffs(p)
     for i in eachindex(cs)
         cs[i] *= l
     end
 end
-#=function _copyto!(x::Vector, y::Vector)
-    resize!(x, length(y))
-    copyto!(x, y)
-    nothing
-end
-function Base.copyto!(x::Vector{MPoly}, y::Vector{MPoly})
-    tsx = terms(x)
-    tsy = terms(y)
-    resize!(tsx, length(tsy))
-    for i in eachindex(tsy)
-        tsx[i] = copy(tsy[i])
-    end
-    x
-end
-function Base.copyto!(x::MPoly, y::MPoly)
-    _copyto!(terms(x), terms(y))
-end
-function Base.copyto!(x::SparsePoly, y::SparsePoly)
-    _copyto!(coeffs(x), coeffs(y))
-    _copyto!(x.exps, y.exps)
-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)
-    # pp = similar(p)
+    p = copy(p)
     while !iszero(p) && degree(p) >= degree(d)
         s = mulpow!(dd, d, lc(p), degree(p) - degree(d))
-        # copyto!(pp, p);
-        # p, pp = pp, p
-        p = copy(p)
         mul!(p, l)
         sub!(p, s)
         k -= 1
@@ -256,32 +236,76 @@ function pseudorem(p::SparsePoly, d::SparsePoly)
     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])
-    # TODO: short circuit and split to small terms
-    for i in 3:length(cfs)
+    isone(g) || for i in 3:length(cfs)
         g = gcd(g, cfs[i])
+        isone(g) && break
     end
     g
 end
 
-function univariate_gcd(x::AbstractPoly, y::AbstractPoly)
-    while !iszero(y)
-        x = pseudorem(x, y)
-        x, y = y, x
+#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 x#, a / x, b / x
+    return q, r
 end
 
-function divexact(a::SparsePoly{T}, b::T) where {T<:AbstractPoly}
-    cfs = MPoly[divexact(c, b) for c in coeffs(a)]
+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, b)
+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
@@ -346,14 +370,14 @@ function emplace_back!(poly, ts, perm, chunk_start_idx, idx, olddegree)
     push!(poly.coeffs, coeff)
     nothing
 end
-term2polycoeff(t::Term, v::IDType) = Term(t.coeff, Monomial(filter(!isequal(v), monomial(t).ids)))
+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
-        count(isequal(v), monomial(t).ids)
+        degree(monomial(t), v)
     end
     perm = sortperm(pows, rev=true)
-    coeffs = MPoly[]
+    coeffs = typeof(p)[]
     exps = UInt[]
     poly = SparsePoly(coeffs, exps, v)
 
@@ -374,7 +398,8 @@ function SparsePoly(p::MPoly, v::IDType)
     return poly
 end
 
-function univariate_to_multivariate(g::SparsePoly{<:AbstractPoly})
+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)
@@ -390,3 +415,17 @@ function univariate_to_multivariate(g::SparsePoly{<:AbstractPoly})
     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/utils.jl b/src/utils.jl
new file mode 100644
index 0000000..c9c1809
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,61 @@
+@generated function reduce_tup(f::F, inds::Tuple{Vararg{Any,N}}) where {F,N}
+    q = Expr(:block, Expr(:meta, :inline, :propagate_inbounds))
+    if N == 1
+        push!(q.args, :(inds[1]))
+        return q
+    end
+    syms = Vector{Symbol}(undef, N)
+    i = 0
+    for n ∈ 1:N
+        syms[n] = iₙ = Symbol(:i_, (i += 1))
+        push!(q.args, Expr(:(=), iₙ, Expr(:ref, :inds, n)))
+    end
+    W =  1 << (8sizeof(N) - 2 - leading_zeros(N))
+    while W > 0
+        _N = length(syms)
+        for _ ∈ 2W:W:_N
+            for w ∈ 1:W
+                new_sym = Symbol(:i_, (i += 1))
+                push!(q.args, Expr(:(=), new_sym, Expr(:call, :f, syms[w], syms[w+W])))
+                syms[w] = new_sym
+            end
+            deleteat!(syms, 1+W:2W)
+        end
+        W >>>= 1
+    end
+    q
+end
+
+struct GetWithShift{F,T}
+    f::F
+    tup::T
+    shift::Int
+end
+(g::GetWithShift)(i) = g.f(g.tup[i + g.shift])
+
+struct GetWithShift2{F,T}
+    f::F
+    tup1::T
+    tup2::T
+end
+(g::GetWithShift2)(i) = g.f(g.tup1[i], g.tup2[i])
+struct GetWithShift3{F,T,S}
+    f::F
+    tup::T
+    s::S
+end
+(g::GetWithShift3)(i) = g.f(g.tup[i], g.s)
+
+function _fmap(f::F, x::Tuple{Vararg{Any,K}}) where {K,F}
+    ntuple(GetWithShift(f, x, 0), Val(K))
+end
+function _fmap(f::F, x::Tuple{Vararg{Any,K}}, y::Tuple{Vararg{Any,K}}) where {K,F}
+    ntuple(GetWithShift2(f, x, y), Val(K))
+end
+function _fmap(f::F, x::Tuple{Vararg{Any,K}}, y) where {K,F}
+    ntuple(GetWithShift3(f, x, y), Val(K))
+end
+
+Base.@pure __parameterless_type(T) = Base.typename(T).wrapper
+parameterless_type(x) = parameterless_type(typeof(x))
+parameterless_type(x::Type) = __parameterless_type(x)
diff --git a/test/runtests.jl b/test/runtests.jl
index baf084b..d1f6e6f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,9 +1,41 @@
 using LoopPoly
+using LoopPoly: divexact
 using LoopPoly: lc
 using Test
 
 LoopPoly.debugmode() = true
 
+@testset "Fuzz SparsePoly" begin
+    x = Uninomial(0, 1)
+    for i in 1:1000
+        p = sum(rand(-2:2)*x^i for i in 0:10)
+        q = sum(rand(-2:2)*x^i for i in 0:5)
+        g = gcd(p, q)
+        if !isone(g)
+            @show g
+        end
+        pdg = divexact(p, g)
+        qdg = divexact(q, g)
+        @test gcd(pdg, qdg) == one(x)
+    end
+end
+
+@testset "Fuzz MPoly" begin
+    x, y, z = [Monomial([i]) for i in 0:2]
+    for _ in 1:1000
+        pterms = [rand(-2:2)*prod(rand([x, y, z])^i for i in 0:3) for j in 0:7]
+        p = sum(pterms)
+        q = sum(rand(-2:2)*prod(rand([x, y, z])^i for i in 0:2)  for j in 0:5)
+        g = gcd(p, q)
+        if !isone(g)
+            @show g
+        end
+        pdg = divexact(p, g)
+        qdg = divexact(q, g)
+        @test gcd(pdg, qdg) == one(x)
+    end
+end
+
 @testset "pseudorem" begin
     x = Uninomial(0, 1)
     p = 2x^10 + x^7 + 7*x^2 + x + 3x
@@ -19,7 +51,7 @@ end
 function test_gcd(x, y)
     g1 = gcd(x, y)
     g2 = gcd(y, x)
-    @test sign(lc(g1)) * g1 == sign(lc(g1)) * g2
+    @test sign(lc(g1)) * g1 == sign(lc(g2)) * g2
     return g1
 end
 
@@ -73,4 +105,34 @@ end
     g = gcd(m, m)
     @test g == m
     @test LoopPoly.degree(g) == 60 + 18 + 12 + 24
+
+    c1 = 10*(x * z + x)
+    c2 = 2*(x^2 + z)
+    c3 = 2*(2 - z  )
+    c4 = 20*(x * z^2)
+    e1 = 0
+    e2 = 5
+    e3 = 7
+    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
+    for i in 0:3
+        @test test_gcd(p + i, q) == p + i
+    end
+
+    k = y^2 + 1
+    @test test_gcd(x*k, z*k) == k
+    @test test_gcd(x*k, (z+1)*k) == k
+    @test test_gcd(x*k, p*k) == k
+
+    p = 2*x*y
+    q = 2*x*y + x
+    @test test_gcd(q, p) == x
+
+    p = x*y + y
+    q = -x*y - y
+    @test test_gcd(q, p) == p
+
+    @test test_gcd((-x + 1) * (y^2+1), -x+1) == -x + 1
 end