Trac #18461: Implement Field._gcd_univariate_polynomial()

The goal of this ticket is to implement a Cython method
`Field._gcd_univariate_polynomial()` using the standard Euclidean
algorithm.  This is much faster than
`Fields.ParentMethods._gcd_univariate_polynomial()`, which calls
`EuclideanDomains.ElementMethods.gcd()`, since both are Python methods.
(The `_gcd_univariate_polynomial` mechanism was introduced in #13442.)

The following bug can then be fixed by just removing
`PolynomialRealDense.gcd()`, which does not take into account the case
where one of the arguments of `gcd` is zero:
{{{
sage: R.<x> = RR[]
sage: x.gcd(R.zero())
Traceback (most recent call last):
...
TypeError: 'MinusInfinity' object cannot be interpreted as an index
}}}
Removing this method ''without'' implementing
`Field._gcd_univariate_polynomial()` (falling back on Python code) is
about twice as slow; with this ticket there is no slowdown.

Similarly, to make `CC` and `QQbar` use the new method instead of Python
code, we also remove the now obsolete `Polynomial_generic_field.gcd()`.

URL: http://trac.sagemath.org/18461
Reported by: pbruin
Ticket author(s): Peter Bruin
Reviewer(s): Bruno Grenet

Author: Release Manager

src/doc/en/thematic_tutorials/coercion_and_categories.rst

@@ -122,6 +122,7 @@ This base class provides a lot more methods than a general parent::
      '_coerce_self',
      '_coerce_try',
      '_default_category',
+     '_gcd_univariate_polynomial',
      '_gens',
      '_has_coerce_map_from',
      '_ideal_class_',

src/sage/rings/polynomial/polynomial_element_generic.py

@@ -726,37 +726,6 @@ def quo_rem(self, other):
             R = R[:R.degree()] - (aaa*B[:B.degree()]).shift(diff_deg)
         return (Q, R)
 
-    @coerce_binop
-    def gcd(self, other):
-        """
-        Return the greatest common divisor of this polynomial and ``other``, as
-        a monic polynomial.
-
-        INPUT:
-
-        - ``other`` -- a polynomial defined over the same ring as ``self``
-
-        EXAMPLES::
-
-            sage: R.<x> = QQbar[]
-            sage: (2*x).gcd(2*x^2)
-            x
-
-            sage: zero = R.zero()
-            sage: zero.gcd(2*x)
-            x
-            sage: (2*x).gcd(zero)
-            x
-            sage: zero.gcd(zero)
-            0
-        """
-        from sage.categories.euclidean_domains import EuclideanDomains
-        g = EuclideanDomains().ElementMethods().gcd(self, other)
-        c = g.leading_coefficient()
-        if c.is_unit():
-            return (1/c)*g
-        return g
-
     @coerce_binop
     def xgcd(self, other):
         r"""

src/sage/rings/polynomial/polynomial_real_mpfr_dense.pyx

@@ -612,49 +612,6 @@ cdef class PolynomialRealDense(Polynomial):
         r._normalize()
         return q, r * leading
 
-    @coerce_binop
-    def gcd(self, other):
-        """
-        Returns the gcd of self and other as a monic polynomial. Due to the
-        inherit instability of division in this inexact ring, the results may
-        not be entirely stable.
-
-        EXAMPLES::
-
-            sage: R.<x> = RR[]
-            sage: (x^3).gcd(x^5+1)
-            1.00000000000000
-            sage: (x^3).gcd(x^5+x^2)
-            x^2
-            sage: f = (x+3)^2 * (x-1)
-            sage: g = (x+3)^5
-            sage: f.gcd(g)
-            x^2 + 6.00000000000000*x + 9.00000000000000
-
-        Unless the division is exact (i.e. no rounding occurs) the returned gcd is
-        almost certain to be 1. ::
-
-            sage: f = (x+RR.pi())^2 * (x-1)
-            sage: g = (x+RR.pi())^5
-            sage: f.gcd(g)
-            1.00000000000000
-
-        """
-        aval = self.valuation()
-        a = self >> aval
-        bval = other.valuation()
-        b = other >> bval
-        if b.degree() > a.degree():
-            a, b = b, a
-        while b: # this will be exactly zero when the previous b is degree 0
-            q, r = a.quo_rem(b)
-            a, b = b, r
-        if a.degree() == 0:
-            # make sure gcd of "relatively prime" things is exactly 1
-            return self._parent(1) << min(aval, bval)
-        else:
-            return a * ~a[a.degree()] << min(aval, bval)
-
     def __call__(self, *args, **kwds):
         """
         EXAMPLES::

src/sage/rings/ring.pyx

@@ -2182,6 +2182,53 @@ cdef class Field(PrincipalIdealDomain):
         """
         raise NotImplementedError, "Algebraic closures of general fields not implemented."
 
+    def _gcd_univariate_polynomial(self, a, b):
+        """
+        Return the gcd of ``a`` and ``b`` as a monic polynomial.
+
+        .. WARNING:
+
+            If the base ring is inexact, the results may not be
+            entirely stable.
+
+        TESTS::
+
+            sage: for A in (RR, CC, QQbar):
+            ....:     g = A._gcd_univariate_polynomial
+            ....:     R.<x> = A[]
+            ....:     z = R.zero()
+            ....:     assert(g(2*x, 2*x^2) == x and
+            ....:            g(z, 2*x) == x and
+            ....:            g(2*x, z) == x and
+            ....:            g(z, z) == z)
+
+            sage: R.<x> = RR[]
+            sage: (x^3).gcd(x^5+1)
+            1.00000000000000
+            sage: (x^3).gcd(x^5+x^2)
+            x^2
+            sage: f = (x+3)^2 * (x-1)
+            sage: g = (x+3)^5
+            sage: f.gcd(g)
+            x^2 + 6.00000000000000*x + 9.00000000000000
+
+        The following example illustrates the fact that for inexact
+        base rings, the returned gcd is often 1 due to rounding::
+
+            sage: f = (x+RR.pi())^2 * (x-1)
+            sage: g = (x+RR.pi())^5
+            sage: f.gcd(g)
+            1.00000000000000
+
+        """
+        while b:
+            q, r = a.quo_rem(b)
+            a, b = b, r
+        if a:
+            a = a.monic()
+        return a
+
+
 cdef class Algebra(Ring):
     """
     Generic algebra