sagemathgh-38924: avoid unnecessary divisions and calls to gcd

    
In an effort to make sagemath#38108 more usable, we optimize the computation of
gcd's in generic polynomial rings.
    
URL: sagemath#38924
Reported by: Martin Rubey
Reviewer(s): Martin Rubey, mathzeta, Travis Scrimshaw

src/sage/arith/misc.py

@@ -2060,7 +2060,7 @@ def xgcd(a, b=None):
         sage: h = R.base_ring().gen()
         sage: S.<y> = R.fraction_field()[]
         sage: xgcd(y^2, a*h*y + b)
-        (1, 7*a^2/b^2, (((-h)*a)/b^2)*y + 1/b)
+        (1, 7*a^2/b^2, (((-7)*a)/(h*b^2))*y + 7/(7*b))
 
     Tests with randomly generated integers::
 

src/sage/categories/fields.py

@@ -673,6 +673,11 @@ def gcd(self, other):
                 sage: gcd(0.0, 0.0)                                                     # needs sage.rings.real_mpfr
                 0.000000000000000
 
+            TESTS::
+
+                sage: QQbar(0).gcd(QQbar.zeta(3))
+                1
+
             AUTHOR:
 
             - Simon King (2011-02) -- :issue:`10771`
@@ -687,7 +692,7 @@ def gcd(self, other):
                 from sage.rings.integer_ring import ZZ
                 try:
                     return P(ZZ(self).gcd(ZZ(other)))
-                except TypeError:
+                except (TypeError, ValueError):
                     pass
 
             if self == P.zero() and other == P.zero():

src/sage/categories/unique_factorization_domains.py

@@ -152,7 +152,7 @@ def _gcd_univariate_polynomial(self, f, g):
                 sage: S.<T> = R[]
                 sage: p = (-3*x^2 - x)*T^3 - 3*x*T^2 + (x^2 - x)*T + 2*x^2 + 3*x - 2
                 sage: q = (-x^2 - 4*x - 5)*T^2 + (6*x^2 + x + 1)*T + 2*x^2 - x
-                sage: quo,rem=p.pseudo_quo_rem(q); quo,rem
+                sage: quo, rem = p.pseudo_quo_rem(q); quo, rem
                 ((3*x^4 + 13*x^3 + 19*x^2 + 5*x)*T + 18*x^4 + 12*x^3 + 16*x^2 + 16*x,
                  (-113*x^6 - 106*x^5 - 133*x^4 - 101*x^3 - 42*x^2 - 41*x)*T
                          - 34*x^6 + 13*x^5 + 54*x^4 + 126*x^3 + 134*x^2 - 5*x - 50)
@@ -167,50 +167,71 @@ def _gcd_univariate_polynomial(self, f, g):
                 sage: g = 4*x + 2
                 sage: f.gcd(g).parent() is R
                 True
+
+            A slightly more exotic base ring::
+
+                sage: P = PolynomialRing(QQbar, 4, "x")
+                sage: p = sum(QQbar.zeta(i + 1) * P.gen(i) for i in range(4))
+                sage: ((p^4 - 1).gcd(p^3 + 1) / (p + 1)).is_unit()
+                True
             """
+            def content(X):
+                """
+                Return the content of ``X`` up to a unit.
+                """
+                # heuristically, polynomials tend to be monic
+                X_it = reversed(X.coefficients())
+                x = next(X_it)
+                # TODO: currently, there is no member of
+                # `UniqueFactorizationDomains` with `is_unit`
+                # significantly slower than `is_one`, see
+                # :issue:`38924` - when such a domain is eventually
+                # implemented, check whether this is a bottleneck
+                if x.is_unit():
+                    return None
+                for c in X_it:
+                    x = x.gcd(c)
+                    if x.is_unit():
+                        return None
+                return x
+
             if f.degree() < g.degree():
-                A,B = g, f
+                A, B = g, f
             else:
-                A,B = f, g
+                A, B = f, g
 
             if B.is_zero():
                 return A
 
-            a = b = self.zero()
-            for c in A.coefficients():
-                a = a.gcd(c)
-                if a.is_one():
-                    break
-            for c in B.coefficients():
-                b = b.gcd(c)
-                if b.is_one():
-                    break
-
-            d = a.gcd(b)
-            A = A // a
-            B = B // b
-            g = h = 1
-
-            delta = A.degree()-B.degree()
-            _,R = A.pseudo_quo_rem(B)
+            a = content(A)
+            if a is not None:
+                A //= a
+            b = content(B)
+            if b is not None:
+                B //= b
 
+            one = A.base_ring().one()
+            if a is not None and b is not None:
+                d = a.gcd(b)
+            else:
+                d = one
+            g = h = one
+            delta = A.degree() - B.degree()
+            _, R = A.pseudo_quo_rem(B)
             while R.degree() > 0:
                 A = B
-                B = R // (g*h**delta)
+                h_delta = h**delta
+                B = R // (g * h_delta)
                 g = A.leading_coefficient()
-                h = h*g**delta // h**delta
+                h = h * g**delta // h_delta
                 delta = A.degree() - B.degree()
                 _, R = A.pseudo_quo_rem(B)
 
             if R.is_zero():
-                b = self.zero()
-                for c in B.coefficients():
-                    b = b.gcd(c)
-                    if b.is_one():
-                        break
-
-                return d*B // b
-
+                b = content(B)
+                if b is None:
+                    return d * B
+                return d * B // b
             return f.parent()(d)
 
     class ElementMethods:

src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py

@@ -613,7 +613,8 @@ def normalize_coordinates(self):
             sage: f.normalize_coordinates(); f
             Dynamical system of Projective Berkovich line over Cp(3) of precision 20
              induced by the map
-              Defn: Defined on coordinates by sending (x : y) to (x^2 : y^2)
+              Defn: Defined on coordinates by sending (x : y)
+               to ((2 + O(3^20))*x^2 : (2 + O(3^20))*y^2)
 
 
         Normalize_coordinates may sometimes fail over `p`-adic fields::

src/sage/schemes/affine/affine_morphism.py

@@ -541,7 +541,7 @@ def homogenize(self, n):
             Scheme endomorphism of Projective Space of dimension 1
              over Algebraic Field
               Defn: Defined on coordinates by sending (x0 : x1) to
-                    (x0*x1 : 1/2*x0^2 + x0*x1 + 3/2*x1^2)
+                    (2*x0*x1 : x0^2 + 2*x0*x1 + 3*x1^2)
 
         ::