Rename _div into div_multiple

Author: blegat

src/division.jl

@@ -44,38 +44,50 @@ _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_multiple(m1::AbstractMonomialLike, m2::AbstractMonomialLike, ::MA.MutableTrait)
+    return mapexponents(-, m1, m2)
+end
+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(t1::AbstractTermLike, t2::AbstractTermLike)
-    term(_div(coefficient(t1), coefficient(t2)), _div(monomial(t1), monomial(t2)))
+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 _div(f::APL, g::APL)
+function div_multiple(f::APL, g::APL, mf=MA.IsNotMutable())
     lt = leadingterm(g)
-    rf = MA.copy_if_mutable(f)
+    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)
     end
     return q
 end
+# 2-arguments useful to be used with `Base.Fix2`
+_div_multiple!(a, b) = div_multiple(a, b, MA.IsMutable())
 
 Base.div(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g; kwargs...)[1]
 Base.rem(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g; kwargs...)[2]
@@ -98,7 +110,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 +148,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 +267,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 +303,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

@@ -397,7 +397,7 @@ _polynomial(ts, state, ::MA.IsNotMutable) = polynomial(ts, state)
 _polynomial(ts, state, ::MA.IsMutable) = polynomial!(ts, state)
 
 """
-    primitive_univariate_gcd!(p::APL, q::APL, algo::AbstractUnivariateGCDAlgorithm)
+    primitive_univariate_gcd!(p::APL, q::APL, algo::AbstractUnivariateGCDAlgorithm, buffer=nothing)
 
 Returns the `gcd` of primitive polynomials `p` and `q` using algorithm `algo`
 which is a subtype of [`AbstractUnivariateGCDAlgorithm`](@ref).
@@ -419,11 +419,19 @@ end
 
 # If `p` and `q` do not have the same type then the local variables `p` and `q`
 # won't be type stable so we create `u` and `v`.
-function primitive_univariate_gcd!(p::APL, q::APL, algo::GeneralizedEuclideanAlgorithm)
+function primitive_univariate_gcd!(
+    p::APL,
+    q::APL,
+    algo::GeneralizedEuclideanAlgorithm,
+    buffer=nothing,
+)
     if maxdegree(p) < maxdegree(q)
-        return primitive_univariate_gcd!(q, p, algo)
+        return primitive_univariate_gcd!(q, p, algo, buffer)
     end
     R = MA.promote_operation(gcd, typeof(p), typeof(q))
+    if isnothing(buffer)
+        buffer = MA.buffer_for(rem_or_pseudo_rem, R, R, typeof(algo))
+    end
     u = convert(R, p)
     v = convert(R, q)
     while true
@@ -436,7 +444,7 @@ function primitive_univariate_gcd!(p::APL, q::APL, algo::GeneralizedEuclideanAlg
         end
 
         d_before = degree(leadingmonomial(u))
-        r = MA.operate!!(rem_or_pseudo_rem, u, v, algo)
+        r = MA.buffered_operate!!(buffer, rem_or_pseudo_rem, u, v, algo)
         d_after = degree(leadingmonomial(r))
         if d_after == d_before
             not_divided_error(u, v)
@@ -485,7 +493,7 @@ function primitive_univariate_gcdx(u0::APL, v0::APL, algo::GeneralizedEuclideanA
 end
 
 
-function primitive_univariate_gcd!(p::APL, q::APL, ::SubresultantAlgorithm)
+function primitive_univariate_gcd!(p::APL, q::APL, ::SubresultantAlgorithm, buffer=nothing)
     error("Not implemented yet")
 end
 
@@ -512,20 +520,35 @@ If the coefficients are not `AbstractFloat`, this
 *Art of computer programming, volume 2: Seminumerical algorithms.*
 Addison-Wesley Professional. Third edition.
 """
-function univariate_gcd(p1::APL{S}, p2::APL{T}, algo::AbstractUnivariateGCDAlgorithm, m1::MA.MutableTrait, m2::MA.MutableTrait) where {S,T}
-    return univariate_gcd(_field_absorb(algebraic_structure(S), algebraic_structure(T)), p1, p2, algo, m1, m2)
-end
-function univariate_gcd(::UFD, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm, m1::MA.MutableTrait, m2::MA.MutableTrait)
+function univariate_gcd(
+    p1::APL{S},
+    p2::APL{T},
+    algo::AbstractUnivariateGCDAlgorithm,
+    m1::MA.MutableTrait,
+    m2::MA.MutableTrait,
+    buffer=nothing
+) where {S,T}
+    return univariate_gcd(_field_absorb(algebraic_structure(S), algebraic_structure(T)), p1, p2, algo, m1, m2, buffer)
+end
+function univariate_gcd(
+    ::UFD,
+    p1::APL,
+    p2::APL,
+    algo::AbstractUnivariateGCDAlgorithm,
+    m1::MA.MutableTrait,
+    m2::MA.MutableTrait,
+    buffer=nothing,
+)
     f1, g1 = primitive_part_content(p1, algo, m1)
     f2, g2 = primitive_part_content(p2, algo, m2)
-    pp = primitive_univariate_gcd!(f1, f2, algo)
+    pp = primitive_univariate_gcd!(f1, f2, algo, buffer)
     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)
+    return mapcoefficients!(Base.Fix2(MA.mul!!, gg), pp, nonzero = true)
 end
 
 function univariate_gcd(::Field, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm, m1::MA.MutableTrait, m2::MA.MutableTrait)
-    return primitive_univariate_gcd!(_copy(p1, m1), _copy(p2, m2), algo)
+    return primitive_univariate_gcd!(_copy(p1, m1), _copy(p2, m2), algo, buffer)
 end
 
 function univariate_gcdx(p1::APL{S}, p2::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {S,T}
@@ -652,5 +675,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 mapcoefficients(Base.Fix2(_div_multiple!, g), p, mutability; nonzero = true), g
 end

src/monomial.jl

@@ -130,7 +130,7 @@ 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))
 

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/allocations.jl

@@ -151,15 +151,22 @@ function _test_univ_gcd(T, algo)
     @polyvar x
     p1 = o * x^2 + 2o * x + o
     p2 = o * x + o
+    buffer = buffer_for(MP.rem_or_pseudo_rem, typeof(p1), typeof(p2), typeof(algo))
     mutable_alloc_test(mutable_copy(p1), 0) do p1
-        MP.univariate_gcd(p1, p2, algo, IsMutable(), IsMutable())
+        MP.primitive_univariate_gcd!(p1, p2, algo, buffer)
+    end
+    if T === Int
+        mutable_alloc_test(mutable_copy(p1), 0) do p1
+            MP.gcd(p1, p2, algo, IsMutable(), IsMutable())
+        end
     end
 end
 
 function test_univ_gcd()
     algo = GeneralizedEuclideanAlgorithm()
-    _test_univ_gcd(Int, algo)
-    #_test_univ_gcd(BigInt, algo) # TODO
+    @testset "$T" for T in [Int, BigInt]
+        _test_univ_gcd(T, algo)
+    end
 end
 
 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