0-allocation mutable rem (#235)

* 0-allocation mutable rem

* Rename constant operations

* Checkout DP#bl/mutable_rem

* Fix constant_function -> right_constant_function

* Add Benchmark 0

* Mutable one in primitive_univariate_gcd

* 0-alloc for BigInt as well

* Format

* Tests fixes

* Fix

* Rename type arg

* Add independence tests

* Fix missing copies

* Benchmark more types

* Exclude test on Julia v1.6

* Missing blank line

* Factor out divided check

* Add test

* Add rem_or_pseudo_rem

* Export

* Fix

* Fix

* Fix

* Mutable rem for Rational{BigInt}

* Fix

* Fix test

* Add TODO

Author: blegat

.github/workflows/ci.yml

@@ -51,7 +51,7 @@ jobs:
         shell: julia --project=@. {0}
         run: |
           using Pkg
-          Pkg.add(PackageSpec(name="DynamicPolynomials", rev="master"))
+          Pkg.add(PackageSpec(name="DynamicPolynomials", rev="bl/mutable_rem"))
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:

docs/src/division.md

@@ -12,4 +12,6 @@ The `divrem` function returns ``(q, r)``.
 
 ```@docs
 divides
+pseudo_rem
+rem_or_pseudo_rem
 ```

perf/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+SIMDPolynomials = "4edc584b-a88b-4acf-80bb-891198a11e01"
+TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"

perf/gcd.jl

@@ -1,30 +1,30 @@
 using BenchmarkTools
 
-function bench0(x, T=Rational{BigInt})
+function bench0(x, T)
     o = one(T)
     p = o * x^2 + 2o * x + o
     q = o * x + o
     @benchmark gcd($p, $q)
 end
 
-function bench1(x, y, z, T=Rational{BigInt})
+function bench1(x, y, z, T)
     o = one(T)
     p = (o * x + o * y^2) * (o * z^3 + o * y^2 + o * x)
     q = (o * x + o * y + o * z) * (o * x^2 + o * y)
     @benchmark gcd($p, $q)
 end
 
-function bench2(x, y, z, t)
-    p = x * y + 3 * (z * t)
-    q = (p + 1) * p
+function bench2(x, y, z, t, T)
+    p = x * y + T(3) * (z * t)
+    q = (p + T(1)) * p
     @benchmark gcd($p, $q)
 end
 
-function bench3(x, y, z, t)
-    c1 = 10*(x * z + x)
-    c2 = 2*(x^2 + z)
-    c3 = 2*(2 - z  )
-    c4 = 20*(x * z^2)
+function bench3(x, y, z, t, T)
+    c1 = T(10)*(x * z + x)
+    c2 = T(2)*(x^2 + z)
+    c3 = T(2)*(2 - z)
+    c4 = T(20)*(x * z^2)
     e1 = 0
     e2 = 5
     e3 = 7
@@ -37,39 +37,63 @@ function bench3(x, y, z, t)
 end
 
 import SIMDPolynomials
-function bench_SP()
+function bench_SP(T)
     x, y, z, t = [SIMDPolynomials.PackedMonomial{4,7}(i) for i in 0:3]
-    b1 = bench1(x, y, z)
-    b2 = bench2(x, y, z, t)
-    b3 = bench3(x, y, z, t)
-    return b1, b2, b3
+    b0 = bench0(x, T)
+    b1 = bench1(x, y, z, T)
+    b2 = bench2(x, y, z, t, T)
+    b3 = bench3(x, y, z, t, T)
+    return b0, b1, b2, b3
 end
 
 import DynamicPolynomials
-function bench_DP()
+function bench_DP(T)
     DynamicPolynomials.@polyvar x y z t
-    b1 = bench1(x, y, z)
-    b2 = bench2(x, y, z, t)
-    b3 = bench3(x, y, z, t)
-    return b1, b2, b3
+    b0 = bench0(x, T)
+    b1 = bench1(x, y, z, T)
+    b2 = bench2(x, y, z, t, T)
+    b3 = bench3(x, y, z, t, T)
+    return b0, b1, b2, b3
 end
 
 import TypedPolynomials
-function bench_TP()
+function bench_TP(T)
     TypedPolynomials.@polyvar x y z t
-    b1 = bench1(x, y, z)
-    b2 = bench2(x, y, z, t)
-    b3 = bench3(x, y, z, t)
-    return b1, b2, b3
+    b0 = bench0(x, T)
+    b1 = bench1(x, y, z, T)
+    b2 = bench2(x, y, z, t, T)
+    b3 = bench3(x, y, z, t, T)
+    return b0, b1, b2, b3
 end
 
 include("table.jl")
 
-function bench()
-    bs1, bs2, bs3 = bench_SP()
-    bd1, bd2, bd3 = bench_DP()
-    bt1, bt2, bt3 = bench_TP()
+function bench(T)
+    bs0, bs1, bs2, bs3 = bench_SP(T)
+    bd0, bd1, bd2, bd3 = bench_DP(T)
+    bt0, bt1, bt2, bt3 = bench_TP(T)
+    println("### Benchmark 0")
+    println()
+    prettyprint(bs0, bd0, bt0)
+    println()
+    println("### Benchmark 1")
+    println()
     prettyprint(bs1, bd1, bt1)
+    println()
+    println("### Benchmark 2")
+    println()
     prettyprint(bs2, bd2, bt2)
+    println()
+    println("### Benchmark 3")
+    println()
     prettyprint(bs3, bd3, bt3)
+    println()
+end
+
+function bench()
+    for T in [Int, BigInt, Rational{BigInt}]
+        println("## $T")
+        println()
+        bench(T)
+    end
 end

src/comparison.jl

@@ -6,7 +6,7 @@ Base.iszero(t::AbstractTerm) = iszero(coefficient(t))
 Base.iszero(t::AbstractPolynomial) = iszero(nterms(t))
 
 # See https://github.com/blegat/MultivariatePolynomials.jl/issues/22
-# avoids the call to be transfered to eqconstant
+# avoids the call to be transfered to left_constant_eq
 Base.:(==)(α::Nothing, x::APL) = false
 Base.:(==)(x::APL, α::Nothing) = false
 Base.:(==)(α::Dict, x::APL) = false
@@ -26,19 +26,19 @@ function polyeqterm(p::AbstractPolynomial, t)
 end
 polyeqterm(p::APL, t) = polyeqterm(polynomial(p), t)
 
-eqconstant(α, v::AbstractVariable) = false
-eqconstant(v::AbstractVariable, α) = false
-function _termeqconstant(t::AbstractTermLike, α)
+left_constant_eq(α, v::AbstractVariable) = false
+right_constant_eq(v::AbstractVariable, α) = false
+function _term_constant_eq(t::AbstractTermLike, α)
     if iszero(t)
         iszero(α)
     else
         α == coefficient(t) && isconstant(t)
     end
 end
-eqconstant(α, t::AbstractTermLike) = _termeqconstant(t, α)
-eqconstant(t::AbstractTermLike, α) = _termeqconstant(t, α)
-eqconstant(α, p::APL) = polyeqterm(p, α)
-eqconstant(p::APL, α) = polyeqterm(p, α)
+left_constant_eq(α, t::AbstractTermLike) = _term_constant_eq(t, α)
+right_constant_eq(t::AbstractTermLike, α) = _term_constant_eq(t, α)
+left_constant_eq(α, p::APL) = polyeqterm(p, α)
+right_constant_eq(p::APL, α) = polyeqterm(p, α)
 
 function Base.:(==)(mono::AbstractMonomial, v::AbstractVariable)
     return isone(degree(mono)) && variable(mono) == v

src/default_polynomial.jl

@@ -150,32 +150,85 @@ end
 Base.copy(p::VectorPolynomial) = MA.mutable_copy(p)
 
 function grlex end
-function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{T,TT}, q::Union{Polynomial,AbstractTermLike}) where {T,TT}
-    get1 = let p=p
-        i -> p.terms[i]
+
+function __polynomial_merge!(op::MA.AddSubMul, p::Polynomial{T,TT}, get1, set, push, resize, keep, q, t::AbstractTermLike...) where {T,TT}
+    compare_monomials = let t=t, get1=get1, q=q
+        (tp, j) -> begin
+            if tp isa Int && j isa Int
+                tp = get1(tp)
+            end
+            grlex(*(monomials(q)[j], monomial.(t)...), monomial(tp))
+        end
     end
-    get2 = let p=p
+    get2 = let t=t, q=q
         i -> begin
-            t = terms(q)[i]
-            TT(MA.scaling_convert(T, MA.operate(op, coefficient(t))), monomial(t))
+            tq = terms(q)[i]
+            TT(MA.scaling_convert(T, MA.operate(MA.add_sub_op(op), *(coefficient(tq), coefficient.(t)...))), *(monomial(tq), monomial.(t)...))
         end
     end
-    set = let p=p
-        (i, t) -> begin
-            p.terms[i] = t
+    combine = let t=t, p=p, q=q
+        (i, j) -> begin
+            if i isa Int
+                p.terms[i] = Term(MA.operate!!(op, coefficient(p.terms[i]), coefficients(q)[j], coefficient.(t)...), monomial(p.terms[i]))
+            else
+                typeof(i)(MA.operate!!(op, coefficient(i), coefficients(q)[j], coefficient.(t)...), monomial(i))
+            end
         end
     end
-    push = let p=p
-        t -> push!(p.terms, t)
+    polynomial_merge!(
+        nterms(p), nterms(q), get1, get2, set, push,
+        compare_monomials, combine, keep, resize
+    )
+    return
+end
+
+function __polynomial_merge!(op::MA.AddSubMul, p::Polynomial{T,TT}, get1, set, push, resize, keep, t::AbstractTermLike, q::Polynomial, buffer=nothing) where {T,TT}
+    compare_monomials = let t=t, get1=get1, q=q
+        (tp, j) -> begin
+            if tp isa Int && j isa Int
+                tp = get1(tp)
+            end
+            grlex(monomial(t) * monomials(q)[j], monomial(tp))
+        end
+    end
+    get2 = let t=t, q=q
+        i -> begin
+            tq = terms(q)[i]
+            TT(MA.scaling_convert(T, MA.operate(MA.add_sub_op(op), coefficient(t) * coefficient(tq))), monomial(t) * monomial(tq))
+        end
+    end
+    combine = let t=t, p=p, q=q
+        (i, j) -> begin
+            if i isa Int
+                p.terms[i] = Term(MA.buffered_operate!!(buffer, op, coefficient(p.terms[i]), coefficient(t), coefficients(q)[j]), monomial(p.terms[i]))
+            else
+                typeof(i)(MA.buffered_operate!!(buffer, op, coefficient(i), coefficient(t), coefficients(q)[j]), monomial(i))
+            end
+        end
     end
-    compare_monomials = let q=q
+    polynomial_merge!(
+        nterms(p), nterms(q), get1, get2, set, push,
+        compare_monomials, combine, keep, resize
+    )
+    return
+end
+
+function __polynomial_merge!(op::Union{typeof(+), typeof(-)}, p::Polynomial{T,TT}, get1, set, push, resize, keep, q::Union{Polynomial,AbstractTermLike}) where {T,TT}
+    compare_monomials = let get1=get1, q=q
         (t, j) -> begin
             if t isa Int && j isa Int
                 t = get1(t)
             end
             grlex(monomials(q)[j], monomial(t))
         end
     end
+    get2 = let q=q
+        i -> begin
+            t = terms(q)[i]
+            # `operate` makes sure we make a copy of the term as it may be stored directly in `p`
+            TT(MA.scaling_convert(T, MA.operate(op, coefficient(t))), monomial(t))
+        end
+    end
     combine = let p=p, q=q
         (i, j) -> begin
             if i isa Int
@@ -185,6 +238,25 @@ function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{T,TT}, q::Un
             end
         end
     end
+    polynomial_merge!(
+        nterms(p), nterms(q), get1, get2, set, push,
+        compare_monomials, combine, keep, resize
+    )
+    return
+end
+
+function _polynomial_merge!(op::Union{typeof(+), typeof(-), MA.AddSubMul}, p::Polynomial{T,TT}, args...) where {T,TT}
+    get1 = let p=p
+        i -> p.terms[i]
+    end
+    set = let p=p
+        (i, t) -> begin
+            p.terms[i] = t
+        end
+    end
+    push = let p=p
+        t -> push!(p.terms, t)
+    end
     resize = let p=p
         (n) -> resize!(p.terms, n)
     end
@@ -198,12 +270,29 @@ function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{T,TT}, q::Un
             end
         end
     end
-    polynomial_merge!(
-        nterms(p), nterms(q), get1, get2, set, push,
-        compare_monomials, combine, keep, resize
-    )
+    __polynomial_merge!(op, p, get1, set, push, resize, keep, args...)
     return p
 end
+
+function MA.operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial, q::Union{AbstractTermLike,Polynomial})
+    return _polynomial_merge!(op, p, q)
+end
+
+function MA.operate!(op::MA.AddSubMul, p::Polynomial, q::Polynomial, args::AbstractTermLike...)
+    return _polynomial_merge!(op, p, q, args...)
+end
+
+function MA.operate!(op::MA.AddSubMul, p::Polynomial, t::AbstractTermLike, q::Polynomial)
+    return _polynomial_merge!(op, p, t, q)
+end
+
+function MA.buffer_for(op::MA.AddSubMul, ::Type{<:Polynomial{S}}, ::Type{<:AbstractTermLike{T}}, ::Type{<:Polynomial{U}}) where {S,T,U}
+    return MA.buffer_for(op, S, T, U)
+end
+function MA.buffered_operate!(buffer, op::MA.AddSubMul, p::Polynomial, t::AbstractTermLike, q::Polynomial)
+    return _polynomial_merge!(op, p, t, q, buffer)
+end
+
 function MA.operate_to!(output::Polynomial, ::typeof(*), p::Polynomial, q::Polynomial)
     empty!(output.terms)
     mul_to_terms!(output.terms, p, q)
@@ -214,6 +303,12 @@ end
 function MA.operate!(::typeof(*), p::Polynomial, q::Polynomial)
     return MA.operate_to!(p, *, MA.mutable_copy(p), q)
 end
+function MA.operate!(::typeof(*), p::Polynomial, t::AbstractTermLike)
+    for i in eachindex(p.terms)
+        p.terms[i] = MA.operate!!(*, p.terms[i], t)
+    end
+    return p
+end
 
 function MA.operate!(::typeof(zero), p::Polynomial)
     empty!(p.terms)
@@ -234,3 +329,9 @@ function MA.operate!(::typeof(removeleadingterm), p::Polynomial)
     pop!(p.terms)
     return p
 end
+
+function MA.operate!(::typeof(unsafe_restore_leading_term), p::Polynomial, t::AbstractTermLike)
+    # We don't need to copy the coefficient of `t`, this is why this function is called `unsafe`
+    push!(p.terms, t)
+    return p
+end

src/default_term.jl

@@ -31,8 +31,8 @@ function MA.promote_operation(::typeof(combine), ::Type{Term{S, M1}}, ::Type{Ter
     return Term{MA.promote_operation(+, S, T), promote_type(M1, M2)}
 end
 
-function MA.mutability(::Type{Term{C,T}}) where {C,T}
-    if MA.mutability(C) isa MA.IsMutable && MA.mutability(T) isa MA.IsMutable
+function MA.mutability(::Type{Term{C,M}}) where {C,M}
+    if MA.mutability(C) isa MA.IsMutable && MA.mutability(M) isa MA.IsMutable
         return MA.IsMutable()
     else
         return MA.IsNotMutable()