Trac #13628: refactor xgcd

Currently, some classes define a _xgcd method which is called by a xgcd
method of euclidean domain elements. Most elements already define xgcd
directly, and I see no real use in having this method in between.

The attached patch renames all _xgcd methods to xgcd and also
streamlines the use of @coerce_binop on them.

This is similar to #13441.

URL: http://trac.sagemath.org/13628
Reported by: saraedum
Ticket author(s): Julian Rueth
Reviewer(s): Niles Johnson, Travis Scrimshaw

src/sage/categories/fields.py

@@ -6,6 +6,7 @@
 #                          William Stein <wstein@math.ucsd.edu>
 #                2008      Teresa Gomez-Diaz (CNRS) <Teresa.Gomez-Diaz@univ-mlv.fr>
 #                2008-2009 Nicolas M. Thiery <nthiery at users.sf.net>
+#                2012      Julian Rueth <julian.rueth@fsfe.org>
 #
 #  Distributed under the terms of the GNU General Public License (GPL)
 #                  http://www.gnu.org/licenses/
@@ -20,6 +21,7 @@
 from sage.misc.cachefunc import cached_method
 from sage.misc.lazy_attribute import lazy_class_attribute
 import sage.rings.ring
+from sage.structure.element import coerce_binop
 
 class Fields(Category_singleton):
     """
@@ -397,3 +399,42 @@ def lcm(self,other):
             if self==0 or other==0:
                 return P.zero()
             return P.one()
+
+        @coerce_binop
+        def xgcd(self, other):
+            """
+            Compute the extended gcd of ``self`` and ``other``.
+
+            INPUT:
+
+            - ``other`` -- an element with the same parent as ``self``
+
+            OUTPUT:
+
+            A tuple ``(r, s, t)`` of elements in the parent of ``self`` such
+            that ``r = s * self + t * other``. Since the computations are done
+            over a field, ``r`` is zero if ``self`` and ``other`` are zero,
+            and one otherwise.
+
+            AUTHORS:
+
+            - Julian Rueth (2012-10-19): moved here from
+              :class:`sage.structure.element.FieldElement`
+
+            EXAMPLES::
+
+                sage: (1/2).xgcd(2)
+                (1, 2, 0)
+                sage: (0/2).xgcd(2)
+                (1, 0, 1/2)
+                sage: (0/2).xgcd(0)
+                (0, 0, 0)
+            """
+            R = self.parent()
+            if not self.is_zero():
+                return (R.one(), ~self, R.zero())
+            if not other.is_zero():
+                return (R.one(), R.zero(), ~other)
+            # else both are 0
+            return (R.zero(), R.zero(), R.zero())
+

src/sage/rings/arith.py

@@ -1884,7 +1884,7 @@ def xgcd(a, b):
         sage: 4*56 + (-5)*44
         4
         sage: g, a, b = xgcd(5/1, 7/1); g, a, b
-        (1, 3, -2)
+        (1, 1/5, 0)
         sage: a*(5/1) + b*(7/1) == g
         True
         sage: x = polygen(QQ)

src/sage/rings/integer.pyx

@@ -181,7 +181,7 @@ import sage.rings.infinity
 import sage.libs.pari.pari_instance
 cdef PariInstance pari = sage.libs.pari.pari_instance.pari
 
-from sage.structure.element import canonical_coercion
+from sage.structure.element import canonical_coercion, coerce_binop
 from sage.misc.superseded import deprecated_function_alias
 
 cdef object numpy_long_interface = {'typestr': '=i4' if sizeof(long) == 4 else '=i8' }
@@ -5304,100 +5304,120 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement):
            return [z, -z]
         return z
 
-    def _xgcd(self, Integer n, bint minimal=0):
+    @coerce_binop
+    def xgcd(self, Integer n):
         r"""
-        Computes extended gcd of self and `n`.
+        Return the extended gcd of this element and ``n``.
 
         INPUT:
 
-        -  ``self`` - integer
+        - ``n`` -- an integer
+
+        OUTPUT:
+
+        A triple ``(g, s, t)`` such that ``g`` is the non-negative gcd of
+        ``self`` and ``n``, and ``s`` and ``t`` are cofactors satisfying the
+        Bezout identity
+
+        .. MATH::
+
+            g = s \cdot \mathrm{self} + t \cdot n.
+
+        .. NOTE::
+
+            There is no guarantee that the cofactors will be minimal. If you
+            need the cofactors to be minimal use :meth:`_xgcd`. Also, using
+            :meth:`_xgcd` directly might be faster in some cases, see
+            :trac:`13628`.
+
+        EXAMPLES::
+
+            sage: 6.xgcd(4)
+            (2, 1, -1)
+
+        """
+        return self._xgcd(n)
+
+    def _xgcd(self, Integer n, bint minimal=0):
+        r"""
+        Return the exteded gcd of ``self`` and ``n``.
 
-        -  ``n`` - integer
+        INPUT:
 
-        -  ``minimal`` - boolean (default false), whether to
-           compute minimal cofactors (see below)
+        - ``n`` -- an integer
+        - ``minimal`` -- a boolean (default: ``False``), whether to compute
+          minimal cofactors (see below)
 
         OUTPUT:
 
-        A triple (g, s, t), where `g` is the non-negative gcd of self
-        and `n`, and `s` and `t` are cofactors satisfying the Bezout identity
+        A triple ``(g, s, t)`` such that ``g`` is the non-negative gcd of
+        ``self`` and ``n``, and ``s`` and ``t`` are cofactors satisfying the
+        Bezout identity
 
-        .. math::
+        .. MATH::
 
-             g = s \cdot \mbox{\rm self} + t \cdot n.
+            g = s \cdot \mathrm{self} + t \cdot n.
 
-        .. note::
+        .. NOTE::
 
-           If minimal is False, then there is no guarantee that the returned
-           cofactors will be minimal in any sense; the only guarantee is that
-           the Bezout identity will be satisfied (see examples below).
+            If ``minimal`` is ``False``, then there is no guarantee that the
+            returned cofactors will be minimal in any sense; the only guarantee
+            is that the Bezout identity will be satisfied (see examples below).
 
-           If minimal is True, the cofactors will satisfy the following
-           conditions. If either self or `n` are zero, the trivial solution
-           is returned. If both self and `n` are nonzero, the function returns
-           the unique solution such that `0 \leq s < |n|/g` (which then must
-           also satisfy `0 \leq |t| \leq |\mbox{\rm self}|/g`).
+            If ``minimal`` is ``True``, the cofactors will satisfy the following
+            conditions. If either ``self`` or ``n`` are zero, the trivial
+            solution is returned. If both ``self`` and ``n`` are nonzero, the
+            function returns the unique solution such that `0 \leq s < |n|/g`
+            (which then must also satisfy
+            `0 \leq |t| \leq |\mbox{\rm self}|/g`).
 
         EXAMPLES::
 
-            sage: Integer(5)._xgcd(7)
+            sage: 5._xgcd(7)
             (1, 3, -2)
             sage: 5*3 + 7*-2
             1
-            sage: g,s,t = Integer(58526524056)._xgcd(101294172798)
+            sage: g,s,t = 58526524056._xgcd(101294172798)
             sage: g
             22544886
             sage: 58526524056 * s + 101294172798 * t
             22544886
 
-        Without minimality guarantees, weird things can happen::
+        Try ``minimal`` option with various edge cases::
 
-            sage: Integer(3)._xgcd(21)
-            (3, 1, 0)
-            sage: Integer(3)._xgcd(24)
-            (3, 1, 0)
-            sage: Integer(3)._xgcd(48)
-            (3, 1, 0)
-
-        Weirdness disappears with minimal option::
-
-            sage: Integer(3)._xgcd(21, minimal=True)
-            (3, 1, 0)
-            sage: Integer(3)._xgcd(24, minimal=True)
-            (3, 1, 0)
-            sage: Integer(3)._xgcd(48, minimal=True)
-            (3, 1, 0)
-            sage: Integer(21)._xgcd(3, minimal=True)
-            (3, 0, 1)
-
-        Try minimal option with various edge cases::
-
-            sage: Integer(5)._xgcd(0, minimal=True)
+            sage: 5._xgcd(0, minimal=True)
             (5, 1, 0)
-            sage: Integer(-5)._xgcd(0, minimal=True)
+            sage: (-5)._xgcd(0, minimal=True)
             (5, -1, 0)
-            sage: Integer(0)._xgcd(5, minimal=True)
+            sage: 0._xgcd(5, minimal=True)
             (5, 0, 1)
-            sage: Integer(0)._xgcd(-5, minimal=True)
+            sage: 0._xgcd(-5, minimal=True)
             (5, 0, -1)
-            sage: Integer(0)._xgcd(0, minimal=True)
+            sage: 0._xgcd(0, minimal=True)
             (0, 1, 0)
 
+        Output may differ with and without the ``minimal`` option::
+
+            sage: 2._xgcd(-2)
+            (2, 1, 0)
+            sage: 2._xgcd(-2, minimal=True)
+            (2, 0, -1)
+
         Exhaustive tests, checking minimality conditions::
 
             sage: for a in srange(-20, 20):
-            ...     for b in srange(-20, 20):
-            ...       if a == 0 or b == 0: continue
-            ...       g, s, t = a._xgcd(b)
-            ...       assert g > 0
-            ...       assert a % g == 0 and b % g == 0
-            ...       assert a*s + b*t == g
-            ...       g, s, t = a._xgcd(b, minimal=True)
-            ...       assert g > 0
-            ...       assert a % g == 0 and b % g == 0
-            ...       assert a*s + b*t == g
-            ...       assert s >= 0 and s < abs(b)/g
-            ...       assert abs(t) <= abs(a)/g
+            ....:   for b in srange(-20, 20):
+            ....:     if a == 0 or b == 0: continue
+            ....:     g, s, t = a._xgcd(b)
+            ....:     assert g > 0
+            ....:     assert a % g == 0 and b % g == 0
+            ....:     assert a*s + b*t == g
+            ....:     g, s, t = a._xgcd(b, minimal=True)
+            ....:     assert g > 0
+            ....:     assert a % g == 0 and b % g == 0
+            ....:     assert a*s + b*t == g
+            ....:     assert s >= 0 and s < abs(b)/g
+            ....:     assert abs(t) <= abs(a)/g
 
         AUTHORS:
 

src/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py

@@ -19,6 +19,7 @@
 import sage.rings.padics.precision_error as precision_error
 import sage.rings.fraction_field_element as fraction_field_element
 import copy
+from sage.structure.element import coerce_binop
 
 from sage.libs.all import pari, pari_gen
 from sage.libs.ntl.all import ZZX
@@ -957,12 +958,52 @@ def _quo_rem_naive(self, right):
     #def lcm(self, right):
     #    raise NotImplementedError
 
+    @coerce_binop
     def xgcd(self, right):
+        """
+        Extended gcd of ``self`` and ``other``.
+
+        INPUT:
+
+        - ``other`` -- an element with the same parent as ``self``
+
+        OUTPUT:
+
+        Polynomials ``g``, ``u``, and ``v`` such that ``g = u*self + v*other``
+
+        .. WARNING::
+
+            The computations are performed using the standard Euclidean
+            algorithm which might produce mathematically incorrect results in
+            some cases. See :trac:`13439`.
+
+        EXAMPLES::
+
+            sage: R.<x> = Qp(3,3)[]
+            sage: f = x + 1
+            sage: f.xgcd(f^2)
+            ((1 + O(3^3))*x + (1 + O(3^3)), (1 + O(3^3)), 0)
+
+        In these examples the results are incorrect, see :trac:`13439`::
+
+            sage: R.<x> = Qp(3,3)[]
+            sage: f = 3*x + 7
+            sage: g = 5*x + 9
+            sage: f.xgcd(f*g) # not tested - currently we get the incorrect result ((O(3^2))*x^2 + (O(3))*x + (1 + O(3^3)), (3^-2 + 2*3^-1 + O(3))*x, (3^-2 + 3^-1 + O(3)))
+            ((3 + O(3^4))*x + (1 + 2*3 + O(3^3)), (1 + O(3^3)), 0)
+
+            sage: R.<x> = Qp(3)[]
+            sage: f = 490473657*x + 257392844/729
+            sage: g = 225227399/59049*x - 8669753175
+            sage: f.xgcd(f*g) # not tested - currently we get the incorrect result ((O(3^18))*x^2 + (O(3^9))*x + (1 + O(3^20)), (3^-5 + 2*3^-1 + 3 + 2*3^2 + 3^4 + 2*3^6 + 3^7 + 3^8 + 2*3^9 + 2*3^11 + 3^13 + O(3^15))*x, (3^5 + 2*3^6 + 2*3^7 + 3^8 + 3^9 + 3^13 + 2*3^14 + 3^17 + 3^18 + 3^20 + 3^21 + 3^23 + 2*3^24 + O(3^25)))
+            ((3^3 + 3^5 + 2*3^6 + 2*3^7 + 3^8 + 2*3^10 + 2*3^11 + 3^12 + 3^13 + 3^15 + 2*3^16 + 3^18 + O(3^23))*x + (2*3^-6 + 2*3^-5 + 3^-3 + 2*3^-2 + 3^-1 + 2*3 + 2*3^2 + 2*3^3 + 2*3^4 + 3^6 + 2*3^7 + 2*3^8 + 2*3^9 + 2*3^10 + 3^11 + O(3^14)), (1 + O(3^20)), 0)
+
+        """
         from sage.misc.stopgap import stopgap
         stopgap("Extended gcd computations over p-adic fields are performed using the standard Euclidean algorithm which might produce mathematically incorrect results in some cases.", 13439)
 
         from sage.rings.polynomial.polynomial_element_generic import Polynomial_generic_field
-        return Polynomial_generic_field._xgcd.im_func(self,right)
+        return Polynomial_generic_field.xgcd(self,right)
 
     #def discriminant(self):
     #    raise NotImplementedError

src/sage/rings/polynomial/polynomial_element_generic.py

@@ -656,7 +656,8 @@ def gcd(self, other):
             return (1/c)*g
         return g
 
-    def _xgcd(self, other):
+    @coerce_binop
+    def xgcd(self, other):
         r"""
         Extended gcd of ``self`` and polynomial ``other``.
 

src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx

@@ -1770,28 +1770,36 @@ cdef class Polynomial_dense_mod_p(Polynomial_dense_mod_n):
         return self.parent()(g, construct=True)
 
     @coerce_binop
-    def xgcd(self, right):
+    def xgcd(self, other):
         r"""
-        Return the extended gcd of self and other, i.e., elements `r, s, t` such that
+        Compute the extended gcd of this element and ``other``.
 
-        .. math::
+        INPUT:
 
-           r = s \cdot self + t \cdot other.
+        - ``other`` -- an element in the same polynomial ring
 
-        """
-        # copied from sage.structure.element.PrincipalIdealDomainElement due to lack of mult inheritance
-        if not PY_TYPE_CHECK(right, Element) or not ((<Element>right)._parent is self._parent):
-            from sage.rings.arith import xgcd
-            return bin_op(self, right, xgcd)
-        return self._xgcd(right)
+        OUTPUT:
 
-    def _xgcd(self, right):
-        """
-        Return ``g, u, v`` such that ``g = u*self + v*right``.
-        """
-        r, s, t = self.ntl_ZZ_pX().xgcd(right.ntl_ZZ_pX())
-        return self.parent()(r, construct=True), self.parent()(s, construct=True), \
-               self.parent()(t, construct=True)
+        A tuple ``(r,s,t)`` of elements in the polynomial ring such
+        that ``r = s*self + t*other``.
+
+        EXAMPLES::
+
+            sage: R.<x> = PolynomialRing(GF(3),implementation='NTL')
+            sage: x.xgcd(x)
+            (x, 0, 1)
+            sage: (x^2 - 1).xgcd(x - 1)
+            (x + 2, 0, 1)
+            sage: R.zero().xgcd(R.one())
+            (1, 0, 1)
+            sage: (x^3 - 1).xgcd((x - 1)^2)
+            (x^2 + x + 1, 0, 1)
+            sage: ((x - 1)*(x + 1)).xgcd(x*(x - 1))
+            (x + 2, 1, 2)
+
+        """
+        r, s, t = self.ntl_ZZ_pX().xgcd(other.ntl_ZZ_pX())
+        return self.parent()(r, construct=True), self.parent()(s, construct=True), self.parent()(t, construct=True)
 
     @coerce_binop
     def resultant(self, other):

src/sage/rings/polynomial/polynomial_quotient_ring.py

@@ -265,7 +265,7 @@ class of the category, and store the current class of the quotient
         sage: first_class == Q.__class__
         False
         sage: [s for s in dir(Q.category().element_class) if not s.startswith('_')]
-        ['cartesian_product', 'gcd', 'is_idempotent', 'is_one', 'is_unit', 'lcm', 'lift']
+        ['cartesian_product', 'gcd', 'is_idempotent', 'is_one', 'is_unit', 'lcm', 'lift', 'xgcd']
 
     As one can see, the elements are now inheriting additional methods: lcm and gcd. Even though
     ``Q.an_element()`` belongs to the old and not to the new element class, it still inherits