Trac #13629: provide xgcd for new polynomial rings through _xgcd_univariate_polynomial

Currently, to add xgcd functionality for a new polynomial ring, one
needs to add a specialized subclass of !PolynomialElement.

The attached patch allows rings to provide a
{{{_xgcd_univariate_polynomial}}} method which will be called by
!PolynomialElement to compute xgcds.

This is similar to #10635, #13442.

URL: http://trac.sagemath.org/13629
Reported by: saraedum
Ticket author(s): Julian Rueth
Reviewer(s): Peter Bruin, Bruno Grenet

Author: Release Manager

src/doc/en/thematic_tutorials/coercion_and_categories.rst

@@ -132,6 +132,7 @@ This base class provides a lot more methods than a general parent::
      '_pseudo_fraction_field',
      '_random_nonzero_element',
      '_unit_ideal',
+     '_xgcd_univariate_polynomial',
      '_zero_element',
      '_zero_ideal',
      'algebraic_closure',

src/sage/categories/fields.py

@@ -207,7 +207,7 @@ def _gcd_univariate_polynomial(self, f, g):
 
             INPUT:
 
-                - ``f``, ``g`` -- two polynomials defined over ``self``
+            - ``f``, ``g`` -- two polynomials defined over ``self``
 
             .. NOTE::
 
@@ -368,6 +368,88 @@ def __pow__(self, n):
             from sage.modules.all import FreeModule
             return FreeModule(self, n)
 
+        def _xgcd_univariate_polynomial(self, left, right):
+            r"""
+            Return an extended gcd of ``left`` and ``right``.
+
+            INPUT:
+
+            - ``left``, ``right`` -- two polynomials over this field
+
+            OUTPUT:
+
+            Polynomials ``g``, ``u``, and ``v`` such that ``g`` is a
+            greatest common divisor of ``left and ``right``, and such
+            that ``g = u*left + v*right`` holds.
+
+            .. NOTE::
+
+                This is a helper method for
+                :meth:`sage.rings.polynomial.polynomial_element.Polynomial.xgcd`.
+
+            EXAMPLES::
+
+                sage: P.<x> = QQ[]
+                sage: F = (x^2 + 2)*x^3; G = (x^2+2)*(x-3)
+                sage: g, u, v = QQ._xgcd_univariate_polynomial(F,G)
+                sage: g, u, v
+                (x^2 + 2, 1/27, -1/27*x^2 - 1/9*x - 1/3)
+                sage: u*F + v*G
+                x^2 + 2
+
+            ::
+
+                sage: g, u, v = QQ._xgcd_univariate_polynomial(x,P(0)); g, u, v
+                (x, 1, 0)
+                sage: g == u*x + v*P(0)
+                True
+                sage: g, u, v = QQ._xgcd_univariate_polynomial(P(0),x); g, u, v
+                (x, 0, 1)
+                sage: g == u*P(0) + v*x
+                True
+
+            TESTS:
+
+            We check that the behavior of xgcd with zero elements is
+            compatible with gcd (:trac:`17671`)::
+
+                sage: R.<x> = QQbar[]
+                sage: zero = R.zero()
+                sage: zero.xgcd(2*x)
+                (x, 0, 1/2)
+                sage: (2*x).xgcd(zero)
+                (x, 1/2, 0)
+                sage: zero.xgcd(zero)
+                (0, 0, 0)
+            """
+            R = left.parent()
+            zero = R.zero()
+            one = R.one()
+            if right.is_zero():
+                if left.is_zero():
+                    return (zero, zero, zero)
+                else:
+                    c = left.leading_coefficient()
+                    return (left/c, one/c, zero)
+            elif left.is_zero():
+                c = right.leading_coefficient()
+                return (right/c, zero, one/c)
+
+            # Algorithm 3.2.2 of Cohen, GTM 138
+            A = left
+            B = right
+            U = one
+            G = A
+            V1 = zero
+            V3 = B
+            while not V3.is_zero():
+                Q, R = G.quo_rem(V3)
+                G, U, V1, V3 = V3, V1, U-V1*Q, R
+            V = (G-A*U)//B
+            lc = G.leading_coefficient()
+            return G/lc, U/lc, V/lc
+
+
     class ElementMethods:
         def euclidean_degree(self):
             r"""

src/sage/rings/polynomial/polynomial_element.pyx

@@ -3774,17 +3774,17 @@ cdef class Polynomial(CommutativeAlgebraElement):
     @coerce_binop
     def gcd(self, other):
         """
-        Compute a greatest common divisor of this polynomial and ``other``.
+        Return a greatest common divisor of this polynomial and ``other``.
 
         INPUT:
 
-            - ``other`` -- a polynomial in the same ring as this polynomial
+        - ``other`` -- a polynomial in the same ring as this polynomial
 
         OUTPUT:
 
-            A greatest common divisor as a polynomial in the same ring as this
-            polynomial. Over a field, the return value will be a monic
-            polynomial.
+        A greatest common divisor as a polynomial in the same ring as
+        this polynomial. If the base ring is a field, the return value
+        is a monic polynomial.
 
         .. NOTE::
 
@@ -6207,6 +6207,66 @@ cdef class Polynomial(CommutativeAlgebraElement):
         """
         return self.parent().variable_name()
 
+    @coerce_binop
+    def xgcd(self, other):
+        r"""
+        Return an extended gcd of this polynomial and ``other``.
+
+        INPUT:
+
+        - ``other`` -- a polynomial in the same ring as this polynomial
+
+        OUTPUT:
+
+        A tuple ``(r, s, t)`` where ``r`` is a greatest common divisor
+        of this polynomial and ``other``, and ``s`` and ``t`` are such
+        that ``r = s*self + t*other`` holds.
+
+        .. NOTE::
+
+            The actual algorithm for computing the extended gcd depends on the
+            base ring underlying the polynomial ring. If the base ring defines
+            a method ``_xgcd_univariate_polynomial``, then this method will be
+            called (see examples below).
+
+        EXAMPLES::
+
+            sage: R.<x> = QQbar[]
+            sage: (2*x^2).gcd(2*x)
+            x
+            sage: R.zero().gcd(0)
+            0
+            sage: (2*x).gcd(0)
+            x
+
+        One can easily add xgcd functionality to new rings by providing a
+        method ``_xgcd_univariate_polynomial``::
+
+            sage: R.<x> = QQ[]
+            sage: S.<y> = R[]
+            sage: h1 = y*x
+            sage: h2 = y^2*x^2
+            sage: h1.xgcd(h2)
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: Univariate Polynomial Ring in x over Rational Field does not provide an xgcd implementation for univariate polynomials
+            sage: T.<x,y> = QQ[]
+            sage: def poor_xgcd(f,g):
+            ....:     ret = S(T(f).gcd(g))
+            ....:     if ret == f: return ret,S.one(),S.zero()
+            ....:     if ret == g: return ret,S.zero(),S.one()
+            ....:     raise NotImplementedError
+            sage: R._xgcd_univariate_polynomial = poor_xgcd
+            sage: h1.xgcd(h2)
+            (x*y, 1, 0)
+            sage: del R._xgcd_univariate_polynomial
+
+        """
+        if hasattr(self.base_ring(), '_xgcd_univariate_polynomial'):
+            return self.base_ring()._xgcd_univariate_polynomial(self, other)
+        else:
+            raise NotImplementedError("%s does not provide an xgcd implementation for univariate polynomials"%self.base_ring())
+
     def variables(self):
         """
         Returns the tuple of variables occurring in this polynomial.

src/sage/rings/polynomial/polynomial_element_generic.py

@@ -726,81 +726,6 @@ def quo_rem(self, other):
             R = R[:R.degree()] - (aaa*B[:B.degree()]).shift(diff_deg)
         return (Q, R)
 
-    @coerce_binop
-    def xgcd(self, other):
-        r"""
-        Extended gcd of ``self`` and polynomial ``other``.
-
-        INPUT:
-
-        - ``other`` -- a polynomial defined over the same ring as ``self``
-
-        OUTPUT:
-
-        Polynomials ``g``, ``u``, and ``v`` such that ``g = u * self + v * other``.
-
-        EXAMPLES::
-
-            sage: P.<x> = QQ[]
-            sage: F = (x^2 + 2)*x^3; G = (x^2+2)*(x-3)
-            sage: g, u, v = F.xgcd(G)
-            sage: g, u, v
-            (x^2 + 2, 1/27, -1/27*x^2 - 1/9*x - 1/3)
-            sage: u*F + v*G
-            x^2 + 2
-
-        ::
-
-            sage: g, u, v = x.xgcd(P(0)); g, u, v
-            (x, 1, 0)
-            sage: g == u*x + v*P(0)
-            True
-            sage: g, u, v = P(0).xgcd(x); g, u, v
-            (x, 0, 1)
-            sage: g == u*P(0) + v*x
-            True
-
-        TESTS:
-
-        We check that the behavior of xgcd with zero elements is compatible with
-        gcd (:trac:`17671`)::
-
-            sage: R.<x> = QQbar[]
-            sage: zero = R.zero()
-            sage: zero.xgcd(2*x)
-            (x, 0, 1/2)
-            sage: (2*x).xgcd(zero)
-            (x, 1/2, 0)
-            sage: zero.xgcd(zero)
-            (0, 0, 0)
-        """
-        R = self.parent()
-        zero = R.zero()
-        one = R.one()
-        if other.is_zero():
-            if self.is_zero():
-                return (zero, zero, zero)
-            else:
-                c = self.leading_coefficient()
-                return (self/c, one/c, zero)
-        elif self.is_zero():
-            c = other.leading_coefficient()
-            return (other/c, zero, one/c)
-
-        # Algorithm 3.2.2 of Cohen, GTM 138
-        A = self
-        B = other
-        U = one
-        G = A
-        V1 = zero
-        V3 = B
-        while not V3.is_zero():
-            Q, R = G.quo_rem(V3)
-            G, U, V1, V3 = V3, V1, U-V1*Q, R
-        V = (G-A*U)//B
-        lc = G.leading_coefficient()
-        return G/lc, U/lc, V/lc
-
 
 class Polynomial_generic_sparse_field(Polynomial_generic_sparse, Polynomial_generic_field):
     """

src/sage/rings/ring.pyx

@@ -2228,6 +2228,68 @@ cdef class Field(PrincipalIdealDomain):
             a = a.monic()
         return a
 
+    def _xgcd_univariate_polynomial(self, a, b):
+        """
+        Return an extended gcd of ``a`` and ``b``.
+
+        INPUT:
+
+        - ``a``, ``b`` -- two univariate polynomials
+
+        OUTPUT:
+
+        A tuple ``(d, u, v)`` of polynomials such that ``d`` is the
+        greatest common divisor (monic or zero) of ``a`` and ``b``,
+        and ``u``, ``v`` satisfy ``d = u*a + v*b``.
+
+        .. WARNING:
+
+            If the base ring is inexact, the results may not be
+            entirely stable.
+
+        ALGORITHM:
+
+        This uses the extended Euclidean algorithm; see for example
+        [Cohen]_, Algorithm 3.2.2.
+
+        REFERENCES:
+
+        .. [Cohen] H. Cohen, A Course in Computational Algebraic
+           Number Theory.  Graduate Texts in Mathematics 138.
+           Springer-Verlag, 1996.
+
+        TESTS::
+
+            sage: for A in (RR, CC, QQbar):
+            ....:     g = A._xgcd_univariate_polynomial
+            ....:     R.<x> = A[]
+            ....:     z, h = R(0), R(1/2)
+            ....:     assert(g(2*x, 2*x^2) == (x, h, z) and
+            ....:            g(z, 2*x) == (x, z, h) and
+            ....:            g(2*x, z) == (x, h, z) and
+            ....:            g(z, z) == (z, z, z))
+
+        """
+        R = a.parent()
+        zero = R.zero()
+        if not b:
+            if not a:
+                return (zero, zero, zero)
+            c = ~a.leading_coefficient()
+            return (c*a, R(c), zero)
+        elif not a:
+            c = ~b.leading_coefficient()
+            return (c*b, zero, R(c))
+        (u, d, v1, v3) = (R.one(), a, zero, b)
+        while v3:
+            q, r = d.quo_rem(v3)
+            (u, d, v1, v3) = (v1, v3, u - v1*q, r)
+        v = (d - a*u) // b
+        if d:
+            c = ~d.leading_coefficient()
+            d, u, v = c*d, c*u, c*v
+        return d, u, v
+
 
 cdef class Algebra(Ring):
     """