avoid unnecessary divisions and calls to gcd

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,49 +167,60 @@ 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
+                True
             """
             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_it = iter(A.coefficients())
+            a = next(A_it)
+            for c in A_it:
                 a = a.gcd(c)
                 if a.is_one():
                     break
-            for c in B.coefficients():
+            else:
+                A //= a
+
+            B_it = iter(B.coefficients())
+            b = next(B_it)
+            for c in B_it:
                 b = b.gcd(c)
                 if b.is_one():
                     break
+            else:
+                B //= b
 
             d = a.gcd(b)
-            A = A // a
-            B = B // b
             g = h = 1
-
-            delta = A.degree()-B.degree()
-            _,R = A.pseudo_quo_rem(B)
-
+            delta = A.degree() - B.degree()
+            _, R = A.pseudo_quo_rem(B)
             while R.degree() > 0:
                 A = B
-                B = R // (g*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_it = iter(B.coefficients())
+                b = next(B_it)
+                for c in B_it:
                     b = b.gcd(c)
                     if b.is_one():
-                        break
-
-                return d*B // b
+                        return d * B
+                return d * B // b  # TODO: does b divide d?
 
             return f.parent()(d)