Trac #16516: Faster roots computation for sparse polynomials over ZZ

Algorithms exist for the computation of the roots of sparse polynomials,
which are much faster than the "generic" algorithms which work for dense
polynomials.

In this ticket, I implement one of these algorithms, for integer roots
of sparse polynomials over ZZ.

URL: https://trac.sagemath.org/16516
Reported by: bruno
Ticket author(s): Bruno Grenet
Reviewer(s): Vincent Delecroix, Travis Scrimshaw, Jeroen Demeyer

Author: Release Manager

src/doc/en/reference/references/index.rst

@@ -396,6 +396,10 @@ REFERENCES:
             the minimum number of vertices", Topology 40 (2001),
             753--772.
 
+.. [CKS1999] Felipe Cucker, Pascal Koiran, and Stephen Smale. *A polynomial-time
+             algorithm for diophantine equations in one variable*, J. Symbolic
+             Computation 27 (1), 21-29, 1999.
+
 .. [CK2015] \J. Campbell and V. Knight. *On testing degeneracy of
             bi-matrix
             games*. http://vknight.org/unpeudemath/code/2015/06/25/on_testing_degeneracy_of_games/ (2015)

src/sage/rings/integer_ring.pyx

@@ -62,6 +62,7 @@ import sage.rings.ideal
 from sage.categories.basic import EuclideanDomains
 from sage.categories.infinite_enumerated_sets import InfiniteEnumeratedSets
 from sage.structure.coerce cimport is_numpy_type
+from sage.structure.element cimport parent
 from sage.structure.parent_gens import ParentWithGens
 from sage.structure.parent cimport Parent
 from sage.structure.richcmp cimport rich_to_bool
@@ -429,9 +430,9 @@ cdef class IntegerRing_class(PrincipalIdealDomain):
         if x in self:
             return self
 
-        from sage.rings.number_field.number_field_element import is_NumberFieldElement
-        if is_NumberFieldElement(x):
-            K, from_K = x.parent().subfield(x)
+        from sage.rings.number_field.number_field_element import NumberFieldElement
+        if isinstance(x, NumberFieldElement):
+            K, from_K = parent(x).subfield(x)
             return K.order(K.gen())
 
         return PrincipalIdealDomain.__getitem__(self, x)
@@ -1216,6 +1217,238 @@ cdef class IntegerRing_class(PrincipalIdealDomain):
         """
         return self(1)
 
+    def _roots_univariate_polynomial(self, p, ring=None, multiplicities=True, algorithm=None):
+        r"""
+        Return the roots of the univariate polynomial ``p``.
+
+        INPUT:
+
+        - ``p`` -- univariate integer polynomial
+
+        - ``ring`` -- ring (default: ``None``); a ring containing `\ZZ` to
+          compute the roots in; ``None`` is equivalent to ``ZZ``
+
+        - ``multiplicities`` -- boolean (default: ``True``); whether to
+          compute the multiplicities
+
+        - ``algorithm`` -- ``"dense"``, ``"sparse"`` or ``None`` (default:
+          ``None``); the algorithm to use
+
+        OUTPUT:
+
+        - If ``multiplicities=True``, the list of pairs `(r, n)` where
+          `r` is a root and `n` the corresponding multiplicity;
+
+        - If ``multiplicities=False``, the list of distincts roots with no
+          information about the multiplicities.
+
+        ALGORITHM:
+
+        If ``algorithm`` is ``"dense"`, the roots are computed using
+        :meth:`_roots_from_factorization`.
+
+        If ``algorithm`` is ``"sparse"``, the roots are computed using the
+        algorithm described in [CKS1999]_.
+
+        If ``algorithm`` is ``None``, use the ``"dense"`` algorithm for
+        polynomials of degree at most `100`, and ``"sparse"`` otherwise.
+
+        .. NOTE::
+
+            This is a helper method for
+            :meth:`sage.rings.polynomial.polynomial_element.Polynomial.roots`.
+
+        TESTS::
+
+            sage: R.<x> = PolynomialRing(ZZ, sparse=True)
+            sage: p = (x + 1)^23 * (x - 1)^23 * (x - 100) * (x + 5445)^5
+            sage: ZZ._roots_univariate_polynomial(p)
+            [(100, 1), (-5445, 5), (1, 23), (-1, 23)]
+            sage: p *= (1 + x^3458645 - 76*x^3435423343 + x^45346567867756556)
+            sage: ZZ._roots_univariate_polynomial(p)
+            [(1, 23), (-1, 23), (100, 1), (-5445, 5)]
+            sage: p *= x^156468451540687043504386074354036574634735074
+            sage: ZZ._roots_univariate_polynomial(p)
+            [(0, 156468451540687043504386074354036574634735074),
+             (1, 23),
+             (-1, 23),
+             (100, 1),
+             (-5445, 5)]
+            sage: ZZ._roots_univariate_polynomial(p, multiplicities=False)
+            [0, 1, -1, 100, -5445]
+
+            sage: R.<x> = PolynomialRing(ZZ, sparse=False)
+            sage: p = (x + 1)^23 * (x - 1)^23 * (x - 100) * (x + 5445)^5
+            sage: ZZ._roots_univariate_polynomial(p)
+            [(100, 1), (-5445, 5), (1, 23), (-1, 23)]
+            sage: ZZ._roots_univariate_polynomial(p, multiplicities=False)
+            [100, -5445, 1, -1]
+
+            sage: ZZ._roots_univariate_polynomial(p, algorithm="sparse")
+            [(100, 1), (-5445, 5), (1, 23), (-1, 23)]
+            sage: ZZ._roots_univariate_polynomial(p, algorithm="dense")
+            [(100, 1), (-5445, 5), (1, 23), (-1, 23)]
+            sage: ZZ._roots_univariate_polynomial(p, algorithm="foobar")
+            Traceback (most recent call last):
+            ...
+            ValueError: unknown algorithm 'foobar'
+
+            sage: p = x^20 * p
+            sage: ZZ._roots_univariate_polynomial(p, algorithm="sparse")
+            [(0, 20), (100, 1), (-5445, 5), (1, 23), (-1, 23)]
+            sage: ZZ._roots_univariate_polynomial(p, algorithm="dense")
+            [(100, 1), (-5445, 5), (0, 20), (1, 23), (-1, 23)]
+        """
+        deg = p.degree()
+        if deg < 0:
+            raise ValueError("roots of 0 are not defined")
+
+        # A specific algorithm is available only for integer roots of integer polynomials
+        if ring is not self and ring is not None:
+            raise NotImplementedError
+
+        # Automatic choice of algorithm
+        if algorithm is None:
+            if deg > 100:
+                algorithm = "sparse"
+            else:
+                algorithm = "dense"
+
+        if algorithm != "dense" and algorithm != "sparse":
+            raise ValueError("unknown algorithm '{}'".format(algorithm))
+
+        # Check if the polynomial is a constant, in which case there are
+        #   no roots. Note that the polynomial is not 0.
+        if deg == 0:
+            return []
+
+        # The dense algorithm is to compute the roots from the factorization.
+        if algorithm == "dense":
+            #NOTE: the content sometimes return an ideal sometimes a number...
+            if parent(p).is_sparse():
+                cont = p.content().gen()
+            else:
+                cont = p.content()
+            if not cont.is_unit():
+                p = p.map_coefficients(lambda c: c // cont)
+            return p._roots_from_factorization(p.factor(), multiplicities)
+
+        v = p.valuation()
+        deg -= v
+        cdef list roots
+
+        # Root 0
+        if v > 0:
+            roots = [(self.zero(), v)] if multiplicities else [self.zero()]
+            if deg == 0: # The shifted polynomial will be constant
+                return roots
+        else:
+            roots = []
+
+        p = p.shift(-v)
+        cdef list e = p.exponents()
+        cdef int i_min, i, j, k = len(e)
+
+        # totally dense polynomial
+        if k == 1 + deg:
+            return roots + p._roots_from_factorization(p.factor(), multiplicities)
+
+        #TODO: the content sometimes return an ideal sometimes a number... this
+        # should be corrected. See
+        # <https://groups.google.com/forum/#!topic/sage-devel/DP_R3rl0vH0>
+        if parent(p).is_sparse():
+            cont = p.content().gen()
+        else:
+            cont = p.content()
+        if not cont.is_unit():
+            p = p.map_coefficients(lambda c: c // cont)
+
+        cdef list c = p.coefficients()
+
+        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+        R = PolynomialRing(p.base_ring(), p.variable_name(), sparse=False)
+        c_max_nbits = c[0].nbits()
+        i_min = 0
+        g = R.zero()
+
+        # Looking for "large" gaps in the exponents
+        # These gaps split the polynomial into lower degree components
+        # Roots of modulus > 1 are common roots of the components
+        for i in range(1, k):
+            if e[i] - e[i-1] > c_max_nbits:
+                g = g.gcd(R( {e[j] - e[i_min]: c[j] for j in range(i_min, i)} ))
+                if g.is_one():
+                    break
+                i_min = i
+                c_max_nbits = c[i].nbits()
+            else:
+                c_max_nbits = max(c[i].nbits(), c_max_nbits)
+
+        # if no gap, directly return the roots of p
+        if g.is_zero():
+            roots.extend(p._roots_from_factorization(p.factor(), multiplicities))
+            return roots
+
+        g = g.gcd(R( {e[j] - e[i_min]: c[j] for j in range(i_min, k)} ))
+
+
+        cdef list cc
+        cdef list ee
+        cdef int m1, m2
+        cdef bint b1, b2
+        # Computation of the roots of modulus 1, without multiplicities
+        # 1 is root iff p(1) == 0 ; -1 iff p(-1) == 0
+        if not multiplicities:
+            if not sum(c):
+                roots.append(self.one())
+            s = 0
+            for j in range(k):
+                s += -c[j] if (e[j] % 2) else c[j]
+            if not s:
+                roots.append(-self.one())
+
+        # Computation of the roots of modulus 1, with multiplicities
+        # For the multiplicities, take the derivatives
+        else:
+            cc = c
+            ee = e
+            m1 = m2 = 0
+            b1 = b2 = True
+
+            for i in range(k):
+                s1 = s2 = 0
+                for j in range(k-i):
+                    if b1: s1 += cc[j]
+                    if b2: s2 += -cc[j] if (ee[j] % 2) else cc[j]
+                if b1 and s1:
+                    m1 = i
+                    b1 = False
+                if b2 and s2:
+                    m2 = i
+                    b2 = False
+                # Stop asap
+                if not (b1 or b2):
+                    break
+
+                # Sparse derivative, that is (p/x^v)' where v = p.val():
+                ee = [ee[j] - ee[0] - 1 for j in range(1,k-i)]
+                cc = [(ee[j] + 1) * cc[j+1] for j in range(k-i-1)]
+
+            if m1 > 0:
+                roots.append((self.one(), m1))
+            if m2 > 0:
+                roots.append((-self.one(), m2))
+
+        # Add roots of modulus > 1 to `roots`:
+        if multiplicities:
+            roots.extend(r for r in g._roots_from_factorization(g.factor(), True)
+                         if r[0].abs() > 1)
+        else:
+            roots.extend(r for r in g._roots_from_factorization(g.factor(), False)
+                         if r.abs() > 1)
+
+        return roots
+
 
     #################################
     ## Coercions to interfaces
@@ -1363,7 +1596,7 @@ def crt_basis(X, xgcd=None):
 
     Y = []
     # 2. Compute extended GCD's
-    ONE = X[0].parent().one()
+    ONE = parent(X[0]).one()
     for i in range(len(X)):
         p = X[i]
         others = P // p

src/sage/rings/polynomial/polynomial_element.pyx

@@ -387,7 +387,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         #         Note that we are guaranteed that right is in the base ring, so this could be fast.
         if not left:
             return self._parent.zero()
-        return self.parent()(left) * self
+        return self._parent(left) * self
 
     cpdef _rmul_(self, Element right):
         """
@@ -407,7 +407,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         #         Note that we are guaranteed that right is in the base ring, so this could be fast.
         if not right:
             return self._parent.zero()
-        return self * self.parent()(right)
+        return self * self._parent(right)
 
     def subs(self, *x, **kwds):
         r"""
@@ -433,7 +433,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             TypeError: keys do not match self's parent
         """
         if len(x) == 1 and isinstance(x[0], dict):
-            g = self.parent().gen()
+            g = self._parent.gen()
             if g in x[0]:
                 return self(x[0][g])
             elif len(x[0]) > 0:
@@ -1279,7 +1279,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.variables()
             (w,)
         """
-        d = dict([(repr(g), R.var(g)) for g in self.parent().gens()])
+        d = dict([(repr(g), R.var(g)) for g in self._parent.gens()])
         return self.subs(**d)
 
     def __invert__(self):
@@ -1293,7 +1293,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: ~f
             1/(x - 90283)
         """
-        return self.parent().one()/self
+        return self.parent().one() / self
 
     def inverse_of_unit(self):
         """
@@ -1304,18 +1304,18 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.inverse_of_unit()
             Traceback (most recent call last):
             ...
-            ValueError: self is not a unit.
+            ValueError: self is not a unit
             sage: f = R(-90283); g = f.inverse_of_unit(); g
             -1/90283
             sage: parent(g)
             Univariate Polynomial Ring in x over Rational Field
         """
         if self.degree() > 0:
             if not self.is_unit():
-                raise ValueError("self is not a unit.")
+                raise ValueError("self is not a unit")
             else:
                 raise NotImplementedError("polynomial inversion over non-integral domains not implemented")
-        return self.parent()(~(self[0]))
+        return self._parent(~(self[0]))
 
     def inverse_mod(a, m):
         """
@@ -1767,7 +1767,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         except NotImplementedError:
             f = self.factor()
 
-        u = self.parent().base_ring()(f.unit())
+        u = self._parent.base_ring()(f.unit())
 
         if all(a[1] % 2 == 0 for a in f) and u.is_square():
             g = u.sqrt()
@@ -1899,7 +1899,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                     raise NotImplementedError
                 allowed_deg_mult = Integer(ring.factored_order()[0][1]) # generally it will be the quotient of this by the degree of the base ring.
             if degree is None:
-                x = self.parent().gen()
+                x = self._parent.gen()
                 if allowed_deg_mult == 1:
                     xq = pow(x,q,self)
                     self = self.gcd(xq-x)
@@ -1954,7 +1954,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                                 break
             if degree == 0:
                 raise ValueError("degree should be nonzero")
-            R = self.parent()
+            R = self._parent
             x = R.gen()
             if degree > 0:
                 xq = x
@@ -2326,12 +2326,12 @@ cdef class Polynomial(CommutativeAlgebraElement):
     def _pow(self, right):
         # TODO: fit __pow__ into the arithmetic structure
         if self.degree() <= 0:
-            return self.parent()(self[0]**right)
+            return self._parent(self[0]**right)
         if right < 0:
             return (~self)**(-right)
-        if (<Polynomial>self) == self.parent().gen():   # special case x**n should be faster!
+        if (<Polynomial>self) == self._parent.gen():   # special case x**n should be faster!
             v = [0]*right + [1]
-            return self.parent()(v, check=True)
+            return self._parent(v, check=True)
         return generic_power(self, right)
 
     def _repr(self, name=None):
@@ -2371,8 +2371,8 @@ cdef class Polynomial(CommutativeAlgebraElement):
         s = " "
         m = self.degree() + 1
         if name is None:
-            name = self.parent().variable_name()
-        atomic_repr = self.parent().base_ring()._repr_option('element_is_atomic')
+            name = self._parent.variable_name()
+        atomic_repr = self._parent.base_ring()._repr_option('element_is_atomic')
         coeffs = self.list(copy=False)
         for n in reversed(xrange(m)):
             x = coeffs[n]
@@ -2458,8 +2458,8 @@ cdef class Polynomial(CommutativeAlgebraElement):
         coeffs = self.list(copy=False)
         m = len(coeffs)
         if name is None:
-            name = self.parent().latex_variable_names()[0]
-        atomic_repr = self.parent().base_ring()._repr_option('element_is_atomic')
+            name = self._parent.latex_variable_names()[0]
+        atomic_repr = self._parent.base_ring()._repr_option('element_is_atomic')
         for n in reversed(xrange(m)):
             x = coeffs[n]
             x = y = latex(x)
@@ -2535,7 +2535,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             {binop:- {binop:** {gen:x {constr_parent: {subscr: {atomic:ZZ}[{atomic:'x'}]} with gens: ('x',)}} {atomic:2}} {atomic:1}}
         """
         if self.degree() > 0:
-            gen = sib.gen(self.parent())
+            gen = sib.gen(self._parent)
             coeffs = self.list(copy=False)
             terms = []
             for i in range(len(coeffs)-1, -1, -1):
@@ -2551,7 +2551,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         elif coerced:
             return sib(self.constant_coefficient(), True)
         else:
-            return sib(self.parent())(sib(self.constant_coefficient(), True))
+            return sib(self._parent)(sib(self.constant_coefficient(), True))
 
     def __setitem__(self, n, value):
         """
@@ -2791,7 +2791,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
 
         - Didier Deshommes (2006-05-25)
         """
-        return self.parent()(polynomial_fateman._mul_fateman_mul(self,right))
+        return self._parent(polynomial_fateman._mul_fateman_mul(self,right))
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -2976,7 +2976,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: (2*x+3).base_ring()
             Integer Ring
         """
-        return self.parent().base_ring()
+        return self._parent.base_ring()
 
     cpdef base_extend(self, R):
         """
@@ -2994,7 +2994,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.change_ring(GF(7))
             x^3 + 4*x + 3
         """
-        S = self.parent().base_extend(R)
+        S = self._parent.base_extend(R)
         return S(self)
 
     def change_variable_name(self, var):
@@ -3010,7 +3010,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.change_variable_name('theta')
             -2/7*theta^3 + 2/3*theta - 19/993
         """
-        R = self.parent().base_ring()[var]
+        R = self._parent.base_ring()[var]
         return R(self.list())
 
     def change_ring(self, R):
@@ -3040,10 +3040,10 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if isinstance(R, Morphism):
             # we're given a hom of the base ring extend to a poly hom
             if R.domain() == self.base_ring():
-                R = self.parent().hom(R, self.parent().change_ring(R.codomain()))
+                R = self._parent.hom(R, self._parent.change_ring(R.codomain()))
             return R(self)
         else:
-            return self.parent().change_ring(R)(self)
+            return self._parent.change_ring(R)(self)
 
     def _mpoly_dict_recursive(self, variables=None, base_ring=None):
         """
@@ -3062,9 +3062,9 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if not self:
             return {}
 
-        var = self.parent().variable_name()
+        var = self._parent.variable_name()
         if variables is None:
-            variables = self.parent().variable_names_recursive()
+            variables = self._parent.variable_names_recursive()
         if not var in variables:
             x = base_ring(self) if base_ring else self
             const_ix = ETuple((0,)*len(variables))
@@ -4039,7 +4039,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                 except NotImplementedError:
                     pass
 
-        R = self.parent().base_ring()
+        R = self._parent.base_ring()
         if hasattr(R, '_factor_univariate_polynomial'):
             return R._factor_univariate_polynomial(self, **kwargs)
 
@@ -4069,7 +4069,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             from_M, to_M = M.structure()
             g = M['x']([to_M(x) for x in self.list()])
             F = g.factor()
-            S = self.parent()
+            S = self._parent
             v = [(S([from_M(x) for x in f.list()]), e) for f, e in F]
             return Factorization(v, from_M(F.unit()))
 
@@ -4107,7 +4107,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                     if R.characteristic() > 1<<29:
                         raise NotImplementedError("Factorization of multivariate polynomials over prime fields with characteristic > 2^29 is not implemented.")
 
-                P = self.parent()
+                P = self._parent
                 P._singular_().set_ring()
                 S = self._singular_().factorize()
                 factors = S[1]
@@ -4166,7 +4166,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         """
         pols = G[0]
         exps = G[1]
-        R = self.parent()
+        R = self._parent
         F = [(R(f), int(e)) for f, e in zip(pols, exps)]
 
         if unit is None:
@@ -4401,7 +4401,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             raise ZeroDivisionError("Pseudo-division by zero is not possible")
 
         # if other is a constant, then R = 0 and Q = self * other^(deg(self))
-        if other in self.parent().base_ring():
+        if other in self._parent.base_ring():
             return (self * other**(self.degree()), self._parent.zero())
 
         R = self
@@ -4413,7 +4413,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         while not R.degree() < B.degree():
             c = R.leading_coefficient()
             diffdeg = R.degree() - B.degree()
-            Q = d*Q + self.parent()(c).shift(diffdeg)
+            Q = d*Q + self._parent(c).shift(diffdeg)
             R = d*R - c*B.shift(diffdeg)
             e -= 1
 
@@ -4669,7 +4669,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             if n is None:
                 q = self.base_ring().order()
                 n = q ** self.degree() - 1
-            y = self.parent().quo(self).gen()
+            y = self._parent.quo(self).gen()
             from sage.groups.generic import order_from_multiple
             return n == order_from_multiple(y, n, n_prime_divs, operation="*")
         else:
@@ -4836,7 +4836,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if sage.rings.rational_field.is_RationalField(R) or is_NumberField(R):
             return NumberField(self, names)
 
-        return R.fraction_field()[self.parent().variable_name()].quotient(self, names)
+        return R.fraction_field()[self._parent.variable_name()].quotient(self, names)
 
     def sylvester_matrix(self, right, variable = None):
         """
@@ -4959,16 +4959,16 @@ cdef class Polynomial(CommutativeAlgebraElement):
         # This code is almost exactly the same as that of
         # sylvester_matrix() in multi_polynomial.pyx.
 
-        if self.parent() != right.parent():
+        if self._parent != right.parent():
             a, b = coercion_model.canonical_coercion(self,right)
             variable = a.parent()(self.variables()[0])
             #We add the variable to cover the case that right is a multivariate
             #polynomial
             return a.sylvester_matrix(b, variable)
 
         if variable:
-            if variable.parent() != self.parent():
-                variable = self.parent()(variable)
+            if variable.parent() != self._parent:
+                variable = self._parent(variable)
 
         from sage.matrix.constructor import matrix
 
@@ -5283,7 +5283,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if self.is_monic():
             return self
         a = ~self.leading_coefficient()
-        R = self.parent()
+        R = self._parent
         if a.parent() != R.base_ring():
             S = R.base_extend(a.parent())
             return a*S(self)
@@ -5306,7 +5306,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.coefficients(sparse=False)
             [1, 0, 2, 0, 1]
         """
-        zero = self.parent().base_ring().zero()
+        zero = self._parent.base_ring().zero()
         if (sparse):
           return [c for c in self.list() if c != zero]
         else:
@@ -5323,7 +5323,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.exponents()
             [0, 2, 4]
         """
-        zero = self.parent().base_ring().zero()
+        zero = self._parent.base_ring().zero()
         l = self.list()
         return [i for i in range(len(l)) if l[i] != zero]
 
@@ -5488,15 +5488,15 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.monomial_coefficient(x^3)
             0
         """
-        if not m.parent() is self.parent():
+        if not m.parent() is self._parent:
             raise TypeError("monomial must have same parent as self.")
 
         d = m.degree()
         coeffs = self.list()
         if 0 <= d < len(coeffs):
             return coeffs[d]
         else:
-            return self.parent().base_ring().zero()
+            return self._parent.base_ring().zero()
 
     def monomials(self):
         """
@@ -5529,8 +5529,8 @@ cdef class Polynomial(CommutativeAlgebraElement):
         """
         if self.is_zero():
             return []
-        v = self.parent().gen()
-        zero = self.parent().base_ring().zero()
+        v = self._parent.gen()
+        zero = self._parent.base_ring().zero()
         coeffs = self.list()
         return [v**i for i in range(self.degree(), -1, -1) if coeffs[i] != zero]
 
@@ -5570,7 +5570,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         """
         n = sage.rings.integer.Integer(n)
         df = self.derivative()
-        K = self.parent().base_ring()
+        K = self._parent.base_ring()
         a = K(x0)
         L = []
         for i in range(n):
@@ -5835,7 +5835,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: pari(a*x + 10*x^3)
             Mod(2, 8)*x^3 + Mod(1, 8)*a*x
         """
-        return self._pari_with_name(self.parent().variable_name())
+        return self._pari_with_name(self._parent.variable_name())
 
     def _pari_or_constant(self, name=None):
         r"""
@@ -5870,7 +5870,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if self.is_constant():
             return self[0].__pari__()
         if name is None:
-            name = self.parent().variable_name()
+            name = self._parent.variable_name()
         return self._pari_with_name(name)
 
     def _pari_with_name(self, name='x'):
@@ -5931,7 +5931,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             a*W^20 + a^7*s*t
         """
         # Get a reference to Magma version of parent.
-        R = magma(self.parent())
+        R = magma(self._parent)
         # Get list of coefficients.
         v = ','.join([a._magma_init_(magma) for a in self.list()])
         return '%s![%s]'%(R.name(), v)
@@ -5977,7 +5977,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.factor()
             (y + 5) * (y + 1)^2
         """
-        R = gap(self.parent())
+        R = gap(self._parent)
         var = list(R.IndeterminatesOfPolynomialRing())[0]
         return self(var)
 
@@ -6102,7 +6102,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         variable = self.variable_name()
         try:
             res = self.__pari__().polresultant(other, variable)
-            return self.parent().base_ring()(res)
+            return self._parent.base_ring()(res)
         except (TypeError, ValueError, PariError, NotImplementedError):
             return self.sylvester_matrix(other).det()
 
@@ -6450,7 +6450,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             (x - c^3) * (x - b^3) * (x - a^3)
 
         """
-        u, v = PolynomialRing(self.parent().base_ring(), ['u', 'v']).gens()
+        u, v = PolynomialRing(self._parent.base_ring(), ['u', 'v']).gens()
         R = (u - v**n).resultant(self(v), v)
         R = R([self.variables()[0], 0])
         if monic:
@@ -6531,7 +6531,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                g = g.monic()
             return g
 
-        fkn = fkd = self.parent().one()
+        fkn = fkd = self._parent.one()
         for j in range(1, k + 1):
             g = star(rpow(self, j), self.symmetric_power(k - j))
             if j % 2:
@@ -6665,7 +6665,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if self.is_zero():
             return self._parent.zero()
         n = self.degree()
-        base_ring = self.parent().base_ring()
+        base_ring = self._parent.base_ring()
         if (is_MPolynomialRing(base_ring) or
             is_PowerSeriesRing(base_ring)):
             # It is often cheaper to compute discriminant of simple
@@ -7325,9 +7325,16 @@ cdef class Polynomial(CommutativeAlgebraElement):
              (0.500000000000000 - 0.866025403784439*I, 5),
              (0.500000000000000 + 0.866025403784439*I, 5)]
         """
-        K = self.parent().base_ring()
+        K = self._parent.base_ring()
+        # If the base ring has a method _roots_univariate_polynomial,
+        # try to use it. An exception is raised if the method does not
+        # handle the current parameters
         if hasattr(K, '_roots_univariate_polynomial'):
-            return K._roots_univariate_polynomial(self, ring=ring, multiplicities=multiplicities, algorithm=algorithm, **kwds)
+            try:
+                return K._roots_univariate_polynomial(self, ring=ring, multiplicities=multiplicities, algorithm=algorithm, **kwds)
+            except NotImplementedError:
+                # This does not handle something, so keep calm and continue on
+                pass
 
         if kwds:
             raise TypeError("roots() got unexpected keyword argument(s): {}".format(kwds.keys()))
@@ -7545,11 +7552,11 @@ cdef class Polynomial(CommutativeAlgebraElement):
             [2, 0, 1/2]
         """
         seq = []
-        K = self.parent().base_ring()
+        K = self._parent.base_ring()
         for fac in F:
             g = fac[0]
             if g.degree() == 1:
-                rt = -g[0]/g[1]
+                rt = -g[0] / g[1]
                 # We need to check that this root is actually in K;
                 # otherwise we'd return roots in the fraction field of K.
                 if rt in K:
@@ -7780,7 +7787,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.variable_name()
             't'
         """
-        return self.parent().variable_name()
+        return self._parent.variable_name()
 
     @coerce_binop
     def xgcd(self, other):
@@ -7924,16 +7931,16 @@ cdef class Polynomial(CommutativeAlgebraElement):
                 if self[k]:
                     return ZZ(k)
         if isinstance(p, Polynomial):
-            p = self.parent().coerce(p)
-        elif is_Ideal(p) and p.ring() is self.parent(): # eventually need to handle fractional ideals in the fraction field
-            if self.parent().base_ring().is_field(): # common case
+            p = self._parent.coerce(p)
+        elif is_Ideal(p) and p.ring() is self._parent: # eventually need to handle fractional ideals in the fraction field
+            if self._parent.base_ring().is_field(): # common case
                 p = p.gen()
             else:
                 raise NotImplementedError
         else:
             from sage.rings.fraction_field import is_FractionField
-            if is_FractionField(p.parent()) and self.parent().has_coerce_map_from(p.parent().ring()):
-                p = self.parent().coerce(p.parent().ring()(p)) # here we require that p be integral.
+            if is_FractionField(p.parent()) and self._parent.has_coerce_map_from(p.parent().ring()):
+                p = self._parent.coerce(p.parent().ring()(p)) # here we require that p be integral.
             else:
                 raise TypeError("The polynomial, p, must have the same parent as self.")
 
@@ -7974,7 +7981,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: f.add_bigoh(2).parent()
             Power Series Ring in x over Integer Ring
         """
-        return self.parent().completion(self.parent().gen())(self).add_bigoh(prec)
+        return self._parent.completion(self._parent.gen())(self).add_bigoh(prec)
 
     @cached_method
     def is_irreducible(self):
@@ -8275,7 +8282,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: del(QQbar._is_squarefree_univariate_polynomial)
 
         """
-        B = self.parent().base_ring()
+        B = self._parent.base_ring()
         if B not in sage.categories.integral_domains.IntegralDomains():
             raise TypeError("is_squarefree() is not defined for polynomials over {}".format(B))
 
@@ -8304,7 +8311,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         # a square-free polynomial has a square-free content
         if not B.is_field():
             content = self.content()
-            if content not in self.parent().base_ring():
+            if content not in self._parent.base_ring():
                 content = content.gen()
             if not content.is_squarefree():
                 return False
@@ -8350,7 +8357,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: (x^3 + x^2).radical()
             x^2 + x
         """
-        P = self.parent()
+        P = self._parent
         R = P.base_ring()
         p = R.characteristic()
         if p == 0 or p > self.degree():
@@ -8542,7 +8549,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             sage: g.parent()
             Univariate Polynomial Ring in x over Finite Field of size 2 (using NTL)
         """
-        R = self.parent()
+        R = self._parent
         if new_base_ring is not None:
             R = R.change_ring(new_base_ring)
         elif isinstance(f, Map):
@@ -8690,7 +8697,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
         if not self.is_monic():
             return False
 
-        P = self.parent()
+        P = self._parent
         gen = P.gen()
 
         if self == gen - 1:  # the first cyc. pol. is treated apart
@@ -8815,7 +8822,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             raise NotImplementedError("not implemented in non-zero characteristic")
         if not S.is_exact():
             raise NotImplementedError("not implemented for inexact base rings")
-        R = self.parent()
+        R = self._parent
         x = R.gen()
         # Extract Phi_n when n is odd.
         t1 = self
@@ -9098,7 +9105,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
             i = self.valuation()
             if i%m:
                 raise ValueError("not a %s power"%m.ordinal_str())
-            S = self.parent()
+            S = self._parent
             return S.gen()**(i//m) * (self>>i).nth_root(m)
         else:
             c = R.characteristic()
@@ -9112,7 +9119,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                         if i%cc:
                             raise ValueError("not a %s power"%m.ordinal_str())
                         ans[i//cc] = self[i].nth_root(cc)
-                p = self.parent()(ans)
+                p = self._parent(ans)
                 m = m // cc
                 if m.is_one():
                     return p
@@ -9174,7 +9181,7 @@ cdef class Polynomial(CommutativeAlgebraElement):
                 raise ValueError("either the dictionary or the specialization must be provided")
         else:
             from sage.rings.polynomial.flatten import SpecializationMorphism
-            phi = SpecializationMorphism(self.parent(),D)
+            phi = SpecializationMorphism(self._parent,D)
         return phi(self)
 
     def _log_series(self, long n):
@@ -10052,13 +10059,13 @@ cdef class Polynomial_generic_dense(Polynomial):
         if self.is_zero():
             return self, self
 
-        R = self.parent().base_ring()
+        R = self._parent.base_ring()
         x = (<Polynomial_generic_dense>self).__coeffs[:] # make a copy
         y = (<Polynomial_generic_dense>other).__coeffs
         m = len(x)  # deg(self)=m-1
         n = len(y)  # deg(other)=n-1
         if m < n:
-            return self.parent().zero(), self
+            return self._parent.zero(), self
 
         quo = list()
         for k from m-n >= k >= 0:

src/sage/rings/polynomial/polynomial_element_generic.py

@@ -141,7 +141,7 @@ def dict(self):
         """
         return dict(self.__coeffs)
 
-    def coefficients(self,sparse=True):
+    def coefficients(self, sparse=True):
         """
         Return the coefficients of the monomials appearing in ``self``.
 
@@ -152,12 +152,22 @@ def coefficients(self,sparse=True):
             7*w^10000 + w^1997 + 5
             sage: f.coefficients()
             [5, 1, 7]
+
+        TESTS:
+
+        Check that all coefficients are in the base ring::
+
+            sage: S.<x> = PolynomialRing(QQ, sparse=True)
+            sage: f = x^4
+            sage: all(c.parent() is QQ for c in f.coefficients(False))
+            True
         """
         if sparse:
-          return [c[1] for c in sorted(six.iteritems(self.__coeffs))]
+            return [self.__coeffs[e] for e in self.exponents()]
         else:
-          return [self.__coeffs[i] if i in self.__coeffs else 0
-                  for i in range(self.degree() + 1)]
+            zero = self.parent().base_ring().zero()
+            return [self.__coeffs[i] if i in self.__coeffs else zero
+                    for i in range(self.degree() + 1)]
 
     def exponents(self):
         """
@@ -190,10 +200,9 @@ def valuation(self):
             sage: R(0).valuation()
             +Infinity
         """
-        c = self.__coeffs.keys()
-        if len(c) == 0:
+        if not self.__coeffs:
             return infinity
-        return ZZ(min(self.__coeffs.keys()))
+        return ZZ(min(self.__coeffs))
 
     def _derivative(self, var=None):
         """
@@ -511,10 +520,9 @@ def degree(self, gen=None):
             sage: f.degree()
             50000
         """
-        v = self.__coeffs.keys()
-        if len(v) == 0:
+        if not self.__coeffs:
             return -1
-        return max(v)
+        return max(self.__coeffs)
 
     def _add_(self, right):
         r"""