Shortcut for dividing by one (#241)

* Specialize on one term

* Export

* Fixes

* DP#bl/mapexponents

* Fixes

* Fix

* Fix

* Fix

Author: blegat

.github/workflows/ci.yml

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

docs/src/division.md

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

src/default_polynomial.jl

@@ -301,14 +301,31 @@ function MA.operate_to!(output::Polynomial, ::typeof(*), p::Polynomial, q::Polyn
     return output
 end
 function MA.operate!(::typeof(*), p::Polynomial, q::Polynomial)
-    return MA.operate_to!(p, *, MA.mutable_copy(p), q)
+    if iszero(q)
+        return MA.operate!(zero, p)
+    elseif nterms(q) == 1
+        return MA.operate!(*, p, leadingterm(q))
+    else
+        return MA.operate_to!(p, *, MA.mutable_copy(p), q)
+    end
 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 mapexponents!(f, p::Polynomial, m::AbstractMonomialLike)
+    for i in eachindex(p.terms)
+        t = p.terms[i]
+        p.terms[i] = term(coefficient(t), mapexponents(f, monomial(t), m))
+    end
+    return p
+end
+function mapexponents(f, p::Polynomial, m::AbstractMonomialLike)
+    P = MA.promote_operation(*, typeof(p), typeof(m))
+    return mapexponents!(f, MA.mutable_copy(convert(P, p)), m)
+end
 
 function MA.operate!(::typeof(zero), p::Polynomial)
     empty!(p.terms)

src/division.jl

@@ -1,4 +1,4 @@
-export divides, pseudo_rem, rem_or_pseudo_rem
+export divides, div_multiple, pseudo_rem, rem_or_pseudo_rem
 
 function Base.round(p::APL; args...)
     # round(0.1) is zero so we cannot use `mapcoefficientsnz`
@@ -44,32 +44,61 @@ _field_absorb(::UFD, ::Field) = Field()
 _field_absorb(::Field, ::UFD) = Field()
 _field_absorb(::Field, ::Field) = Field()
 
-# _div(a, b) assumes that b divides a
-_div(::Field, a, b) = a / b
-_div(::UFD, a, b) = div(a, b)
-_div(a, b) = _div(algebraic_structure(promote_type(typeof(a), typeof(b))), a, b)
-_div(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(-, m1, m2)
-function _div(t::AbstractTerm, m::AbstractMonomial)
-    term(coefficient(t), _div(monomial(t), m))
+"""
+    div_multiple(a, b, ma::MA.MutableTrait)
+
+Return the division of `a` by `b` assuming that `a` is a multiple of `b`.
+If `a` is not a multiple of `b` then this function may return anything.
+"""
+div_multiple(::Field, a, b, ma::MA.MutableTrait) = a / b
+div_multiple(::UFD, a, b, ma::MA.IsMutable) = MA.operate!!(div, a, b)
+div_multiple(::UFD, a, b, ma::MA.IsNotMutable) = div(a, b)
+function div_multiple(a, b, ma::MA.MutableTrait=MA.IsNotMutable())
+    return div_multiple(algebraic_structure(promote_type(typeof(a), typeof(b))), a, b, ma)
 end
-function _div(t1::AbstractTermLike, t2::AbstractTermLike)
-    term(_div(coefficient(t1), coefficient(t2)), _div(monomial(t1), monomial(t2)))
+function div_multiple(m1::AbstractMonomialLike, m2::AbstractMonomialLike, ::MA.MutableTrait=MA.IsNotMutable())
+    return mapexponents(-, m1, m2)
 end
-function _div(f::APL, g::APL)
+function div_multiple(t::AbstractTerm, m::AbstractMonomial, mt::MA.MutableTrait=MA.IsNotMutable())
+    term(_copy(coefficient(t), mt), div_multiple(monomial(t), m))
+end
+function div_multiple(t1::AbstractTermLike, t2::AbstractTermLike, m1::MA.MutableTrait=MA.IsNotMutable())
+    term(div_multiple(coefficient(t1), coefficient(t2), m1), div_multiple(monomial(t1), monomial(t2)))
+end
+function right_constant_div_multiple(f::APL, g, mf::MA.MutableTrait=MA.IsNotMutable())
+    if isone(g)
+        return _copy(f, mf)
+    end
+    return mapcoefficients(coef -> div_multiple(coef, g, mf), f, mf; nonzero = true)
+end
+function div_multiple(f::APL, g::AbstractMonomialLike, mf::MA.MutableTrait=MA.IsNotMutable())
+    if isconstant(g)
+        return _copy(f, mf)
+    end
+    return mapexponents(-, f, g, mf)
+end
+function div_multiple(f::APL, g::AbstractTermLike, mf::MA.MutableTrait=MA.IsNotMutable())
+    f = right_constant_div_multiple(f, coefficient(g), mf)
+    return div_multiple(f, monomial(g), MA.IsMutable())
+end
+function div_multiple(f::APL, g::APL, mf::MA.MutableTrait=MA.IsNotMutable())
     lt = leadingterm(g)
-    rf = MA.copy_if_mutable(f)
+    if nterms(g) == 1
+        return div_multiple(f, lt, mf)
+    end
+    rf = _copy(f, mf)
     rg = removeleadingterm(g)
     q = zero(rf)
     while !iszero(rf)
         ltf = leadingterm(rf)
         if !divides(lt, ltf)
             # In floating point arithmetics, it may happen
             # that `rf` is not zero even if it cannot be reduced further.
-            # As `_div` assumes that `g` divides `f`, we know that
+            # As `div_multiple` assumes that `g` divides `f`, we know that
             # `rf` is approximately zero anyway.
             break
         end
-        qt = _div(ltf, lt)
+        qt = div_multiple(ltf, lt)
         q = MA.add!!(q, qt)
         rf = MA.operate!!(removeleadingterm, rf)
         rf = MA.operate!!(MA.sub_mul, rf, qt, rg)
@@ -98,7 +127,7 @@ function _pseudo_divrem(::UFD, f::APL, g::APL, algo)
     else
         st = constantterm(coefficient(ltg), f)
         new_f = st * removeleadingterm(f)
-        qt = term(coefficient(ltf), _div(monomial(ltf), monomial(ltg)))
+        qt = term(coefficient(ltf), div_multiple(monomial(ltf), monomial(ltg)))
         new_g = qt * rg
         # Check with `::` that we don't have any type unstability on this variable.
         return convert(typeof(f), st), convert(typeof(f), qt), (new_f - new_g)::typeof(f)
@@ -136,11 +165,11 @@ end
 
 function _prepare_s_poly!(::typeof(pseudo_rem), f, ltf, ltg)
     MA.operate!(right_constant_mult, f, coefficient(ltg))
-    return term(coefficient(ltf), _div(monomial(ltf), monomial(ltg)))
+    return term(coefficient(ltf), div_multiple(monomial(ltf), monomial(ltg)))
 end
 
 function _prepare_s_poly!(::typeof(rem), ::APL, ltf, ltg)
-    return _div(ltf, ltg)
+    return div_multiple(ltf, ltg)
 end
 
 function MA.operate!(op::Union{typeof(rem), typeof(pseudo_rem)}, f::APL, g::APL, algo)
@@ -255,7 +284,7 @@ function Base.divrem(f::APL{T}, g::APL{S}; kwargs...) where {T, S}
         if isapproxzero(ltf; kwargs...)
             rf = MA.operate!!(removeleadingterm, rf)
         elseif divides(lm, ltf)
-            qt = _div(ltf, lt)
+            qt = div_multiple(ltf, lt)
             q = MA.add!!(q, qt)
             rf = MA.operate!!(removeleadingterm, rf)
             rf = MA.operate!!(MA.sub_mul, rf, qt, rg)
@@ -291,7 +320,7 @@ function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T
         divisionoccured = false
         for i in useful
             if divides(lm[i], ltf)
-                qt = _div(ltf, lt[i])
+                qt = div_multiple(ltf, lt[i])
                 q[i] = MA.add!!(q[i], qt)
                 rf = MA.operate!!(removeleadingterm, rf)
                 rf = MA.operate!!(MA.sub_mul, rf, qt, rg[i])

src/gcd.jl

@@ -390,7 +390,7 @@ function flatten_variable!(::Type{TT}, poly::APL) where {TT<:AbstractTerm}
             push!(ts, _t * m)
         end
     end
-    return polynomial!(ts)
+    return polynomial!(ts, UniqState())
 end
 
 _polynomial(ts, state, ::MA.IsNotMutable) = polynomial(ts, state)
@@ -521,7 +521,10 @@ function univariate_gcd(::UFD, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgo
     pp = primitive_univariate_gcd!(f1, f2, algo)
     gg = _gcd(g1, g2, algo, MA.IsMutable(), MA.IsMutable())#::MA.promote_operation(gcd, typeof(g1), typeof(g2))
     # Multiply each coefficient by the gcd of the contents.
-    return mapcoefficients!(Base.Fix1(*, gg), pp, nonzero = true)
+    if !isone(gg)
+        MA.operate!(right_constant_mult, pp, gg)
+    end
+    return pp
 end
 
 function univariate_gcd(::Field, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm, m1::MA.MutableTrait, m2::MA.MutableTrait)
@@ -652,5 +655,5 @@ Addison-Wesley Professional. Third edition.
 """
 function primitive_part_content(p, algo::AbstractUnivariateGCDAlgorithm, mutability::MA.MutableTrait)
     g = content(p, algo, MA.IsNotMutable())
-    return mapcoefficients(Base.Fix2(_div, g), p, mutability; nonzero = true), g
+    return right_constant_div_multiple(p, g, mutability), g
 end

src/monomial.jl

@@ -130,13 +130,16 @@ If ``m_1 = \\prod x^{\\alpha_i}`` and ``m_2 = \\prod x^{\\beta_i}`` then it retu
 
 ### Examples
 
-The multiplication `m1 * m2` is equivalent to `mapexponents(+, m1, m2)`, the unsafe division `_div(m1, m2)` is equivalent to `mapexponents(-, m1, m2)`, `gcd(m1, m2)` is equivalent to `mapexponents(min, m1, m2)`, `lcm(m1, m2)` is equivalent to `mapexponents(max, m1, m2)`.
+The multiplication `m1 * m2` is equivalent to `mapexponents(+, m1, m2)`, the unsafe division `div_multiple(m1, m2)` is equivalent to `mapexponents(-, m1, m2)`, `gcd(m1, m2)` is equivalent to `mapexponents(min, m1, m2)`, `lcm(m1, m2)` is equivalent to `mapexponents(max, m1, m2)`.
 """
 mapexponents(f, m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(f, monomial(m1), monomial(m2))
 
 function mapexponents_to! end
 function mapexponents! end
 
+mapexponents(f, a, b, ::MA.IsMutable) = mapexponents!(f, a, b)
+mapexponents(f, a, b, ::MA.IsNotMutable) = mapexponents(f, a, b)
+
 Base.one(::Type{TT}) where {TT<:AbstractMonomialLike} = constantmonomial(TT)
 Base.one(t::AbstractMonomialLike) = constantmonomial(t)
 MA.promote_operation(::typeof(one), MT::Type{<:AbstractMonomialLike}) = monomialtype(MT)

src/operators.jl

@@ -271,6 +271,8 @@ Base.:*(m::AbstractMonomialLike, t::AbstractTermLike) = term(coefficient(t), m *
 Base.:*(t::AbstractTermLike, m::AbstractMonomialLike) = term(coefficient(t), monomial(t) * m)
 Base.:*(t1::AbstractTermLike, t2::AbstractTermLike) = term(coefficient(t1) * coefficient(t2), monomial(t1) * monomial(t2))
 
+MA.operate!(::typeof(*), p::APL, t::AbstractMonomialLike) = mapexponents!(+, p, t)
+Base.:*(p::APL, t::AbstractMonomialLike) = mapexponents(+, p, t)
 Base.:*(t::AbstractTermLike, p::APL) = polynomial!(map(te -> t * te, terms(p)))
 Base.:*(p::APL, t::AbstractTermLike) = polynomial!(map(te -> te * t, terms(p)))
 Base.:*(p::APL, q::APL) = polynomial(p) * polynomial(q)
@@ -361,7 +363,7 @@ function MA.buffered_operate_to!(buffer::AbstractPolynomial, output::AbstractPol
     product = MA.operate_to!!(buffer, *, y, z, args...)
     return MA.operate_to!(output, MA.add_sub_op(op), x, product)
 end
-function MA.buffered_operate!(buffer::AbstractPolynomial, op::MA.AddSubMul, x::AbstractPolynomial, y, z, args::Vararg{Any, N}) where N
+function MA.buffered_operate!(buffer, op::MA.AddSubMul, x::AbstractPolynomial, y, z, args::Vararg{Any, N}) where N
     product = MA.operate_to!!(buffer, *, y, z, args...)
     return MA.operate!(MA.add_sub_op(op), x, product)
 end

src/term.jl

@@ -102,7 +102,7 @@ function coefficient(f::APL, m::AbstractMonomialLike, vars)
         end
         match || continue
 
-        coeff += term(coefficient(t), _div(monomial(t), m))
+        coeff += term(coefficient(t), div_multiple(monomial(t), m))
     end
     coeff
 end

test/monomial.jl

@@ -81,7 +81,7 @@ end
 
     @test monic(x^2) == x^2
 
-    @test MP._div(2x^2*y[1]^3, x*y[1]^2) == 2x*y[1]
+    @test MP.div_multiple(2x^2*y[1]^3, x*y[1]^2) == 2x*y[1]
 
     @test transpose(x) == x
     @test adjoint(x) == x