Implement univariate extended gcd (#178)

* Implement univariate extended gcd

* Fix pseudo division with f becoming zero

* Fixes

Author: blegat

src/division.jl

@@ -56,6 +56,31 @@ end
 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]
 
+function pseudo_divrem(f::APL{S}, g::APL{T}, algo) where {S,T}
+    return _pseudo_divrem(algebraic_structure(S, T), f, g, algo)
+end
+
+function _pseudo_divrem(::Field, f::APL, g::APL, algo)
+    q, r = divrem(f, g)
+    return one(q), q, r
+end
+
+function _pseudo_divrem(::UFD, f::APL, g::APL, algo)
+    ltg = leadingterm(g)
+    rg = removeleadingterm(g)
+    ltf = leadingterm(f)
+    if iszero(f) || !divides(monomial(ltg), ltf)
+        return one(f), zero(f), zero(f)
+    else
+        st = constantterm(coefficient(ltg), f)
+        new_f = st * removeleadingterm(f)
+        qt = term(coefficient(ltf), _div(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)
+    end
+end
+
 function pseudo_rem(f::APL{S}, g::APL{T}, algo) where {S,T}
     return _pseudo_rem(algebraic_structure(S, T), f, g, algo)
 end
@@ -71,7 +96,7 @@ function _pseudo_rem(::UFD, f::APL, g::APL, algo)
     if !divides(monomial(ltg), ltf)
         return false, f
     end
-    while divides(monomial(ltg), ltf)
+    while !iszero(f) && divides(monomial(ltg), ltf)
         new_f = constantterm(coefficient(ltg), f) * removeleadingterm(f)
         new_g = term(coefficient(ltf), _div(monomial(ltf), monomial(ltg))) * rg
         # Check with `::` that we don't have any type unstability on this variable.

src/gcd.jl

@@ -237,6 +237,35 @@ function deflated_gcd(p1::APL{T}, p2::APL{S}, algo) where {T, S}
     end
 end
 
+function Base.gcdx(p1::APL{T}, p2::APL{S}, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) where {T, S}
+    i1, i2, num_common = _extracted_variable(p1, p2)
+    R = MA.promote_operation(gcd, typeof(p1), typeof(p2))
+    if iszero(i1)
+        if iszero(i2)
+            return univariate_gcdx(p1, p2, algo)
+        else
+            if isapproxzero(p1)
+                return zero(R), one(R), convert(R, p2)
+            end
+            error("Not implemented yet")
+        end
+    else
+        if iszero(i2)
+            if isapproxzero(p2)
+                return one(R), zero(R), convert(R, p1)
+            end
+            error("Not implemented yet")
+        else
+            if num_common > 1
+                @assert i1 == i2
+                error("Not implemented yet")
+            else
+                return univariate_gcdx(p1, p2, algo)
+            end
+        end
+    end
+end
+
 # Returns first element in the union of two decreasing vectors
 function _extracted_variable(p1, p2)
     v1 = variables(p1)
@@ -324,6 +353,18 @@ which is a subtype of [`AbstractUnivariateGCDAlgorithm`](@ref).
 """
 function primitive_univariate_gcd end
 
+function not_divided_error(u, v)
+    error(
+        "Polynomial `$v` of degree `$(maxdegree(v))` and effective",
+        " variables `$(effective_variables(v))` does not divide",
+        " polynomial `$u` of degree `$(maxdegree(u))` and effective",
+        " variables `$(effective_variables(u))`. Did you call",
+        " `univariate_gcd` with polynomials with more than one",
+        " variable in common ? If yes, call `gcd` instead, otherwise,",
+        " please report this.",
+    )
+end
+
 # If `p` and `q` do not have the same time then the local variables `p` and `q`
 # won't be type stable.
 function primitive_univariate_gcd(p::APL, q::APL, algo::GeneralizedEuclideanAlgorithm)
@@ -343,15 +384,7 @@ function primitive_univariate_gcd(p::APL, q::APL, algo::GeneralizedEuclideanAlgo
         end
         divided, r = pseudo_rem(u, v, algo)
         if !divided
-            error(
-                "Polynomial `$v` of degree `$(maxdegree(v))` and effective",
-                " variables `$(effective_variables(v))` does not divide",
-                " polynomial `$u` of degree `$(maxdegree(u))` and effective",
-                " variables `$(effective_variables(u))`. Did you call",
-                " `univariate_gcd` with polynomials with more than one",
-                " variable in common ? If yes, call `gcd` instead, otherwise,",
-                " please report this.",
-            )
+            not_divided_error(u, v)
         end
         if !algo.primitive_rem
             r = primitive_part(r, algo)::R
@@ -360,6 +393,42 @@ function primitive_univariate_gcd(p::APL, q::APL, algo::GeneralizedEuclideanAlgo
     end
 end
 
+function primitive_univariate_gcdx(u0::APL, v0::APL, algo::GeneralizedEuclideanAlgorithm)
+    if maxdegree(u0) < maxdegree(v0)
+        a, b, g = primitive_univariate_gcdx(v0, u0, algo)
+        return b, a, g
+    end
+    R = MA.promote_operation(gcd, typeof(u0), typeof(v0))
+    u = convert(R, u0)
+    v = convert(R, v0)
+    if isapproxzero(v)
+        return one(R), zero(R), u
+    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 zero(R), one(R), v
+    end
+    # p * u = q * v + r
+    p, q, r = pseudo_divrem(u, v, algo)
+    if iszero(q)
+        not_divided_error(u, v)
+    end
+    if iszero(r)
+        # Shortcut, does not change the output
+        return zero(R), one(R), v
+    end
+    # TODO
+    #if !algo.primitive_rem
+    #    r = primitive_part(r, algo)::R
+    #end
+    # a * v + b * r = g
+    # a * v + b * (p * u - q * v) = g
+    # b * p * u + (a - b * q) * v = g
+    a, b, g = primitive_univariate_gcdx(v, r, algo)
+    return p * b, (a - b * q), g
+end
+
+
 function primitive_univariate_gcd(p::APL, q::APL, ::SubresultantAlgorithm)
     error("Not implemented yet")
 end
@@ -403,6 +472,21 @@ function univariate_gcd(::Field, p1::APL, p2::APL, algo::AbstractUnivariateGCDAl
     return primitive_univariate_gcd(p1, p2, algo)
 end
 
+function univariate_gcdx(p1::APL{S}, p2::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {S,T}
+    return univariate_gcdx(algebraic_structure(S, T), p1, p2, algo)
+end
+function univariate_gcdx(::UFD, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm)
+    f1, g1 = primitive_part_content(p1, algo)
+    f2, g2 = primitive_part_content(p2, algo)
+    a, b, pp = primitive_univariate_gcdx(f1, f2, algo)
+    gg = _gcd(g1, g2, algo)#::MA.promote_operation(gcd, typeof(g1), typeof(g2))
+    # Multiply each coefficient by the gcd of the contents.
+    return g2 * a, g1 * b, g1 * g2 * mapcoefficientsnz(Base.Fix1(*, gg), pp)
+end
+function univariate_gcdx(::Field, p1::APL, p2::APL, algo::AbstractUnivariateGCDAlgorithm)
+    return primitive_univariate_gcdx(p1, p2, algo)
+end
+
 _gcd(a::APL, b::APL, algo) = gcd(a, b, algo)
 _gcd(a, b::APL, algo) = gcd(a, b, algo)
 _gcd(a::APL, b, algo) = gcd(a, b, algo)

test/division.jl

@@ -96,15 +96,33 @@ function test_gcd_unit(expected, p1, p2, algo)
     @test g == expected || g == -expected
 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)[2])
+    @test a * p1 + b * p2 == g
+end
+function _test_gcdx_unit(expected, p1, p2, algo)
+    test_gcdx_unit(expected, p1, p2, algo)
+    test_gcdx_unit(expected, p2, p1, algo)
+end
+
 
 function univariate_gcd_test(algo=GeneralizedEuclideanAlgorithm())
     Mod.@polyvar x
-    test_gcd_unit(x + 1, x^2 - 1, x^2 + 2x + 1, algo)
-    test_gcd_unit(x + 1, x^2 + 2x + 1, x^2 - 1, algo)
-    test_gcd_unit(x + 1, x^2 - 1, x + 1, algo)
-    test_gcd_unit(x + 1, x + 1, x^2 - 1, algo)
-    test_gcd_unit(x + 1, x - x, x + 1, algo)
-    test_gcd_unit(x + 1, x + 1, x - x, algo)
+    test_gcdx_unit(x + 1, x^2 - 1, x^2 + 2x + 1, algo)
+    test_gcdx_unit(x + 1, x^2 + 2x + 1, x^2 - 1, algo)
+    test_gcdx_unit(x + 1, x^2 - 1, x + 1, algo)
+    test_gcdx_unit(x + 1, x + 1, x^2 - 1, algo)
+    test_gcdx_unit(x + 1, x - x, x + 1, algo)
+    test_gcdx_unit(x + 1, x + 1, x - x, algo)
+    test_gcdx_unit(x - x + 1, x + 1, x + 2, algo)
+    test_gcdx_unit(x - x + 1, x + 1, 2x + 1, algo)
+    test_gcdx_unit(x - x + 1, x - 1, 2x^2 - 2x - 2, algo)
     @test       0  == @inferred gcd(x - x, x^2 - x^2, algo)
     @test       0  == @inferred gcd(x^2 - x^2, x - x, algo)
 end