Add Subresultant algorithm (#264)

* Add Subresultant algorithm

* Fix format

* Exclude subres for gcd

Author: blegat

src/division.jl

@@ -231,6 +231,13 @@ function MA.operate!(
     return MA.buffered_operate!(nothing, op, f, g, algo)
 end
 
+# TODO As suggested in [Knu14, Algorithm R, p. 426] (univariate case only), if
+#      `deg(f) = n` and `deg(g) = m`, during the first `n - m - t` steps, the
+#      coefficient of the `t`th power will simply be multiplied by the leading
+#      coefficient of `g` so we could speed up by multiplying this coefficient
+#      to the power `n - m - t` using `Base.power_by_squaring`.
+#      However, since there might be missing terms, so we don't know in advance
+#      the needed power but we could keep track of it.
 function MA.buffered_operate!(
     buffer,
     op::Union{typeof(rem),typeof(pseudo_rem)},
@@ -240,8 +247,9 @@ function MA.buffered_operate!(
 )
     ltg = leading_term(g)
     ltf = leading_term(f)
+    # This only makes sense in the univariate case but it's only used for univariate gcd anyway
+    skipped_divisions = maxdegree(f) - maxdegree(g) + 1
     MA.operate!(remove_leading_term, g)
-    not_divided_terms = nothing
     while !iszero(f)
         if isapproxzero(ltf) # TODO `, kwargs...)`
             MA.operate!(remove_leading_term, f)
@@ -253,26 +261,24 @@ function MA.buffered_operate!(
             if monomial(ltg) > monomial(ltf)
                 break
             end
-            if isnothing(not_divided_terms)
-                not_divided_terms = term_type(f)
-            end
-            push!(not_divided_terms, ltf)
             MA.operate!(remove_leading_term, f)
         else
             MA.operate!(remove_leading_term, f)
             t = _prepare_s_poly!(op, f, ltf, ltg)
+            skipped_divisions -= 1
             MA.buffered_operate!(buffer, MA.sub_mul, f, t, g)
         end
-        if op === pseudo_rem && algo.primitive_rem
+        if op === pseudo_rem && _primitive_rem(algo)
             f = primitive_part(f, algo, MA.IsMutable())::typeof(f)
         end
-        if algo.skip_last && maxdegree(f) == maxdegree(g)
+        if _skip_last(algo) && maxdegree(f) == maxdegree(g)
             break
         end
         ltf = leading_term(f)
     end
     # Add it back as we cannot modify `g`
     MA.operate!(unsafe_restore_leading_term, g, ltg)
+    _set_skipped_divisions!(algo, skipped_divisions)
     return f
 end
 

src/gcd.jl

@@ -20,6 +20,11 @@ Algorithm computing the greatest common divisor of univariate polynomials using
 the Euclidean algorithm generalized for polynomials with coefficients over a
 a unique factorization domain, see [Knu14, Algorithm E, p. 426-427].
 
+If `primitive_rem` is `true`, the intermediate remainders produced in the
+polynomial division are made primitive. If `primitive_part` is set to `false`,
+only the resuting remainder is made primitive (the intermediate remainders
+of the generalized Euclidean algorithm still need to be made primitive).
+
 [Knu14] Knuth, D.E., 2014.
 *Art of computer programming, volume 2: Seminumerical algorithms.*
 Addison-Wesley Professional. Third edition.
@@ -35,17 +40,41 @@ struct GeneralizedEuclideanAlgorithm <: AbstractUnivariateGCDAlgorithm
     end
 end
 
+_primitive_rem(algo::GeneralizedEuclideanAlgorithm) = algo.primitive_rem
+_skip_last(algo::GeneralizedEuclideanAlgorithm) = algo.skip_last
+function _set_skipped_divisions!(::GeneralizedEuclideanAlgorithm, ::Int) end
+
 """
-    struct SubresultantAlgorithm <: AbstractUnivariateGCDAlgorithm end
+    mutable struct SubresultantAlgorithm <: AbstractUnivariateGCDAlgorithm
+        skipped_divisions::Int
+    end
 
 Algorithm computing the greatest common divisor of univariate polynomials using
 the Subresultant algorithm, see [Knu14, Algorithm C, p. 428-429].
 
+The division by `g*h^δ` in the algorithm only works if the iteration of
+[Knu14, Algorithm R, p. 426] is carried out even when the divided polynomial
+has a zero term. For computational savings, we don't do that so we store
+in `skipped_division` the number of skipped divisions so that the division
+by `g*h^δ` can be adapted accordingly.
+
+In [Knu14, Algorithm C, p. 426], it is stated that there should be ``
+
 [Knu14] Knuth, D.E., 2014.
 *Art of computer programming, volume 2: Seminumerical algorithms.*
 Addison-Wesley Professional. Third edition.
 """
-struct SubresultantAlgorithm <: AbstractUnivariateGCDAlgorithm end
+mutable struct SubresultantAlgorithm <: AbstractUnivariateGCDAlgorithm
+    skipped_divisions::Int
+    SubresultantAlgorithm() = new(0)
+end
+
+_primitive_rem(::SubresultantAlgorithm) = false
+_skip_last(::SubresultantAlgorithm) = false
+function _set_skipped_divisions!(algo::SubresultantAlgorithm, n::Int)
+    algo.skipped_divisions = n
+    return
+end
 
 _coefficient_gcd(α, β) = gcd(α, β)
 _coefficient_gcd(α::AbstractFloat, β) = one(Base.promote_typeof(α, β))
@@ -57,23 +86,23 @@ end
 function Base.lcm(
     p::_APL,
     q::_APL,
-    algo::AbstractUnivariateGCDAlgorithm = GeneralizedEuclideanAlgorithm(),
+    algo::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
 )
     return p * div(q, gcd(p, q, algo))
 end
 function Base.gcd(
     α,
     p::_APL,
-    algo::AbstractUnivariateGCDAlgorithm = GeneralizedEuclideanAlgorithm(),
-    mα::MA.MutableTrait = MA.IsNotMutable(),
+    algo::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
+    ::MA.MutableTrait = MA.IsNotMutable(),
     mp::MA.MutableTrait = MA.IsNotMutable(),
 )
     return _coefficient_gcd(α, content(p, algo, mp))
 end
 function Base.gcd(
     p::_APL,
     α,
-    algo::AbstractUnivariateGCDAlgorithm = GeneralizedEuclideanAlgorithm(),
+    algo::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
     mp::MA.MutableTrait = MA.IsNotMutable(),
     mα::MA.MutableTrait = MA.IsNotMutable(),
 )
@@ -87,7 +116,7 @@ function MA.promote_operation(
     ::typeof(gcd),
     P::Type{<:_APL},
     Q::Type{<:_APL},
-    A::Type = GeneralizedEuclideanAlgorithm,
+    A::Type = SubresultantAlgorithm,
 )
     return MA.promote_operation(rem_or_pseudo_rem, P, Q, A)
 end
@@ -140,7 +169,7 @@ Addison-Wesley Professional. Third edition.
 function Base.gcd(
     p1::_APL{T},
     p2::_APL{S},
-    algo::AbstractUnivariateGCDAlgorithm = GeneralizedEuclideanAlgorithm(),
+    algo::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
     m1::MA.MutableTrait = MA.IsNotMutable(),
     m2::MA.MutableTrait = MA.IsNotMutable(),
 ) where {T,S}
@@ -179,7 +208,7 @@ end
 function Base.gcd(
     t1::AbstractTermLike{T},
     t2::AbstractTermLike{S},
-    algo::AbstractUnivariateGCDAlgorithm = GeneralizedEuclideanAlgorithm(),
+    ::AbstractUnivariateGCDAlgorithm = SubresultantAlgorithm(),
     m1::MA.MutableTrait = MA.IsNotMutable(),
     m2::MA.MutableTrait = MA.IsNotMutable(),
 ) where {T,S}
@@ -558,6 +587,98 @@ function primitive_univariate_gcd!(
     end
 end
 
+function _pow_no_copy(a, b)
+    if isone(b)
+        # `a^1` is `copy(a)` but that copy is not needed here
+        return a
+    else
+        return a^b
+    end
+end
+
+function primitive_univariate_gcd!(
+    p::_APL,
+    q::_APL,
+    algo::SubresultantAlgorithm,
+)
+    if maxdegree(p) < maxdegree(q)
+        return primitive_univariate_gcd!(q, p, algo)
+    end
+    R = MA.promote_operation(gcd, typeof(p), typeof(q))
+    u = convert(R, p)
+    v = convert(R, q)
+    if isapproxzero(v)
+        return primitive_part(u, algo, MA.IsMutable())::R
+    elseif isconstant(v)
+        # `p` and `q` are primitive so if one of them is constant, it cannot
+        # divide the content of the other one.
+        return MA.operate!(one, u)
+    end
+    g = h = nothing # `nothing` means `1`
+    while true
+        δ = maxdegree(u) - maxdegree(v)
+        d_before = degree(leading_monomial(u))
+        r = MA.operate!(rem_or_pseudo_rem, u, v, algo)
+        if isapproxzero(r)
+            return primitive_part(v, algo, MA.IsMutable())::R
+        elseif isconstant(r)
+            return MA.operate!(one, v)
+        end
+
+        d_after = degree(leading_monomial(r))
+        if d_after == d_before
+            not_divided_error(u, v)
+        end
+        if !isnothing(g)
+            if isnothing(h) # equivalent to `iszero(δ)`
+                ghδ = g
+            else
+                ghδ = g * _pow_no_copy(h, δ)
+            end
+            if !iszero(algo.skipped_divisions)
+                @assert algo.skipped_divisions > 0
+                if isnothing(h)
+                    # It is an alias to `g`
+                    ghδ = MA.copy_if_mutable(ghδ)
+                end
+                # TODO not sure this works, sometimes it `ghδ` is not a multiple of `leading_coefficient(v)`
+                #      we just know it divides the content of `r`, it is not guaranteed to be equal to the content of `r`
+                #      We could maybe do better than multiply `r` here though but let's start with this approach as a baseline
+                #ghδ = div_multiple(ghδ, _pow_no_copy(leading_coefficient(v), algo.skipped_divisions), MA.IsMutable())
+                r = MA.operate!(
+                    right_constant_mult,
+                    r,
+                    _pow_no_copy(
+                        leading_coefficient(v),
+                        algo.skipped_divisions,
+                    ),
+                )
+            end
+            r = right_constant_div_multiple(r, ghδ, MA.IsMutable())::R
+        end
+        u, v = v, r::R
+        g = leading_coefficient(u)
+        # Computes `h = h^(1 - δ) * g^δ` (step C3) of [Knu14, Algorithm C p. 429]
+        # If `δ` is zero then `h^(1 - δ) * g^δ = h` so there is nothing to do
+        if δ == 1
+            h = g
+        elseif δ > 1
+            if isnothing(h) || δ == 2
+                # `h^1` is `copy(h)` but that copy is not needed here
+                hδ = h
+            else
+                hδ = h^(δ - 1)
+            end
+            if isnothing(h)
+                h = g
+            else
+                # We assume that `g^δ` is mutable since `δ > 1`
+                h = div_multiple(g^δ, hδ, MA.IsMutable())
+            end
+        end
+    end
+end
+
 function primitive_univariate_gcdx(
     u0::_APL,
     v0::_APL,
@@ -597,7 +718,7 @@ function primitive_univariate_gcdx(
     return p * b, (a - b * q), g
 end
 
-function primitive_univariate_gcd!(p::_APL, q::_APL, ::SubresultantAlgorithm)
+function primitive_univariate_gcdx(p::_APL, q::_APL, ::SubresultantAlgorithm)
     return error("Not implemented yet")
 end
 

test/division.jl

@@ -1,8 +1,7 @@
 using LinearAlgebra, Test
 using Combinatorics
 
-import MutableArithmetics
-const MA = MutableArithmetics
+import MutableArithmetics as MA
 
 using MultivariatePolynomials
 const MP = MultivariatePolynomials
@@ -123,13 +122,15 @@ end
 
 function test_gcdx_unit(expected, p1, p2, algo)
     test_gcd_unit(expected, p1, p2, algo)
-    a, b, g = @inferred gcdx(p1, p2, algo)
-    # it does not make sense, in general, to speak of "the" greatest common
-    # divisor of u and v; there is a set of greatest common divisors, each
-    # one being a unit multiple of the others [Knu14, p. 424] so `expected` and
-    # `-expected` are both accepted.
-    @test iszero(MP.pseudo_rem(g, expected, algo))
-    @test a * p1 + b * p2 == g
+    if !(algo isa SubresultantAlgorithm) # FIXME not implemented yet
+        a, b, g = @inferred gcdx(p1, p2, algo)
+        # it does not make sense, in general, to speak of "the" greatest common
+        # divisor of u and v; there is a set of greatest common divisors, each
+        # one being a unit multiple of the others [Knu14, p. 424] so `expected` and
+        # `-expected` are both accepted.
+        @test iszero(MP.pseudo_rem(g, expected, algo))
+        @test a * p1 + b * p2 == g
+    end
 end
 function _test_gcdx_unit(expected, p1, p2, algo)
     test_gcdx_unit(expected, p1, p2, algo)
@@ -214,7 +215,7 @@ end
 
 function multivariate_gcd_test(
     ::Type{T},
-    algo = GeneralizedEuclideanAlgorithm(),
+    algo = SubresultantAlgorithm(),
 ) where {T}
     Mod.@polyvar x y z
     o = one(T)
@@ -411,6 +412,7 @@ end
     @testset "Division by multiple polynomials examples" begin
         multi_div_test()
     end
+    univariate_gcd_test(SubresultantAlgorithm())
     @testset "Univariate gcd primitive_rem=$primitive_rem" for primitive_rem in
                                                                [false, true]
         @testset "skip_last=$skip_last" for skip_last in [false, true]
@@ -422,6 +424,7 @@ end
     @testset "Multivariate gcd $T" for T in
                                        [Int, BigInt, Rational{BigInt}, Float64]
         if T != Rational{BigInt} || VERSION >= v"1.6"
+            multivariate_gcd_test(T, SubresultantAlgorithm())
             # `gcd` for `Rational{BigInt}` got defined at some point between v1.0 and v1.6
             @testset "primitive_rem=$primitive_rem" for primitive_rem in
                                                         [false, true]