def roots_cyclotomic(f, factor=False):
    """Compute roots of cyclotomic polynomials. """
    L, U = _inv_totient_estimate(f.degree())

    for n in xrange(L, U + 1):
        g = cyclotomic_poly(n, f.gen, polys=True)

        if f == g:
            break
    else:  # pragma: no cover
        raise RuntimeError("failed to find index of a cyclotomic polynomial")

    roots = []

    if not factor:
        for k in xrange(1, n + 1):
            if igcd(k, n) == 1:
                roots.append(exp(2*k*S.Pi*I/n).expand(complex=True))
    else:
        g = Poly(f, extension=(-1)**Rational(1, n))

        for h, _ in g.factor_list()[1]:
            roots.append(-h.TC())

    return sorted(roots, key=default_sort_key)

def roots_cyclotomic(f, factor=False):
    """Compute roots of cyclotomic polynomials. """
    L, U = _inv_totient_estimate(f.degree())

    for n in range(L, U + 1):
        g = cyclotomic_poly(n, f.gen, polys=True)

        if f == g:
            break
    else:  # pragma: no cover
        raise RuntimeError("failed to find index of a cyclotomic polynomial")

    roots = []

    if not factor:
        # get the indices in the right order so the computed
        # roots will be sorted
        h = n//2
        ks = [i for i in range(1, n + 1) if igcd(i, n) == 1]
        ks.sort(key=lambda x: (x, -1) if x <= h else (abs(x - n), 1))
        d = 2*I*pi/n
        for k in reversed(ks):
            roots.append(exp(k*d).expand(complex=True))
    else:
        g = Poly(f, extension=root(-1, n))

        for h, _ in ordered(g.factor_list()[1]):
            roots.append(-h.TC())

    return roots

def random_poly(x, n, inf, sup, domain=ZZ, polys=False):
    """Return a polynomial of degree ``n`` with coefficients in ``[inf, sup]``. """
    poly = Poly(dup_random(n, inf, sup, domain), x, domain=domain)

    if not polys:
        return poly.as_expr()
    else:
        return poly

 def as_poly(self, x=None):
     """Create a Poly instance from ``self``. """
     if x is not None:
         return Poly.new(self.rep, x)
     else:
         if self.alias is not None:
             return Poly.new(self.rep, self.alias)
         else:
             return PurePoly.new(self.rep, Dummy('x'))

def _minimal_polynomial_sq(p, n, x):
    """
    Returns the minimal polynomial for the ``nth-root`` of a sum of surds
    or ``None`` if it fails.

    Parameters
    ==========

    p : sum of surds
    n : positive integer
    x : variable of the returned polynomial

    Examples
    ========

    >>> from sympy.polys.numberfields import _minimal_polynomial_sq
    >>> from sympy import sqrt
    >>> from sympy.abc import x
    >>> q = 1 + sqrt(2) + sqrt(3)
    >>> _minimal_polynomial_sq(q, 3, x)
    x**12 - 4*x**9 - 4*x**6 + 16*x**3 - 8

    """
    from sympy.simplify.simplify import _is_sum_surds

    p = sympify(p)
    n = sympify(n)
    r = _is_sum_surds(p)
    if not n.is_Integer or not n > 0 or not _is_sum_surds(p):
        return None
    pn = p**Rational(1, n)
    # eliminate the square roots
    p -= x
    while 1:
        p1 = _separate_sq(p)
        if p1 is p:
            p = p1.subs({x:x**n})
            break
        else:
            p = p1

    # _separate_sq eliminates field extensions in a minimal way, so that
    # if n = 1 then `p = constant*(minimal_polynomial(p))`
    # if n > 1 it contains the minimal polynomial as a factor.
    if n == 1:
        p1 = Poly(p)
        if p.coeff(x**p1.degree(x)) < 0:
            p = -p
        p = p.primitive()[1]
        return p
    # by construction `p` has root `pn`
    # the minimal polynomial is the factor vanishing in x = pn
    factors = factor_list(p)[1]

    result = _choose_factor(factors, x, pn)
    return result

def _eval_sum_hyper(f, i, a):
    """ Returns (res, cond). Sums from a to oo. """
    from sympy.functions import hyper
    from sympy.simplify import hyperexpand, hypersimp, fraction, simplify
    from sympy.polys.polytools import Poly, factor
    from sympy.core.numbers import Float

    if a != 0:
        return _eval_sum_hyper(f.subs(i, i + a), i, 0)

    if f.subs(i, 0) == 0:
        if simplify(f.subs(i, Dummy('i', integer=True, positive=True))) == 0:
            return S(0), True
        return _eval_sum_hyper(f.subs(i, i + 1), i, 0)

    hs = hypersimp(f, i)
    if hs is None:
        return None

    if isinstance(hs, Float):
        from sympy.simplify.simplify import nsimplify
        hs = nsimplify(hs)

    numer, denom = fraction(factor(hs))
    top, topl = numer.as_coeff_mul(i)
    bot, botl = denom.as_coeff_mul(i)
    ab = [top, bot]
    factors = [topl, botl]
    params = [[], []]
    for k in range(2):
        for fac in factors[k]:
            mul = 1
            if fac.is_Pow:
                mul = fac.exp
                fac = fac.base
                if not mul.is_Integer:
                    return None
            p = Poly(fac, i)
            if p.degree() != 1:
                return None
            m, n = p.all_coeffs()
            ab[k] *= m**mul
            params[k] += [n/m]*mul

    # Add "1" to numerator parameters, to account for implicit n! in
    # hypergeometric series.
    ap = params[0] + [1]
    bq = params[1]
    x = ab[0]/ab[1]
    h = hyper(ap, bq, x)

    return f.subs(i, 0)*hyperexpand(h), h.convergence_statement

def horner(f, *gens, **args):
    """
    Rewrite a polynomial in Horner form.

    Among other applications, evaluation of a polynomial at a point is optimal
    when it is applied using the Horner scheme ([1]).

    Examples
    ========

    >>> from sympy.polys.polyfuncs import horner
    >>> from sympy.abc import x, y, a, b, c, d, e

    >>> horner(9*x**4 + 8*x**3 + 7*x**2 + 6*x + 5)
    x*(x*(x*(9*x + 8) + 7) + 6) + 5

    >>> horner(a*x**4 + b*x**3 + c*x**2 + d*x + e)
    e + x*(d + x*(c + x*(a*x + b)))

    >>> f = 4*x**2*y**2 + 2*x**2*y + 2*x*y**2 + x*y

    >>> horner(f, wrt=x)
    x*(x*y*(4*y + 2) + y*(2*y + 1))

    >>> horner(f, wrt=y)
    y*(x*y*(4*x + 2) + x*(2*x + 1))

    References
    ==========
    [1] - http://en.wikipedia.org/wiki/Horner_scheme

    """
    allowed_flags(args, [])

    try:
        F, opt = poly_from_expr(f, *gens, **args)
    except PolificationFailed as exc:
        return exc.expr

    form, gen = S.Zero, F.gen

    if F.is_univariate:
        for coeff in F.all_coeffs():
            form = form*gen + coeff
    else:
        F, gens = Poly(F, gen), gens[1:]

        for coeff in F.all_coeffs():
            form = form*gen + horner(coeff, *gens, **args)

    return form

def _minpoly_pow(ex, pw, x, dom, mp=None):
    """
    Returns ``minpoly(ex**pw, x)``

    Parameters
    ==========

    ex : algebraic element
    pw : rational number
    x : indeterminate of the polynomial
    dom: ground domain
    mp : minimal polynomial of ``p``

    Examples
    ========

    >>> from sympy import sqrt, QQ, Rational
    >>> from sympy.polys.numberfields import _minpoly_pow, minpoly
    >>> from sympy.abc import x, y
    >>> p = sqrt(1 + sqrt(2))
    >>> _minpoly_pow(p, 2, x, QQ)
    x**2 - 2*x - 1
    >>> minpoly(p**2, x)
    x**2 - 2*x - 1
    >>> _minpoly_pow(y, Rational(1, 3), x, QQ.frac_field(y))
    x**3 - y
    >>> minpoly(y**Rational(1, 3), x)
    x**3 - y

    """
    pw = sympify(pw)
    if not mp:
        mp = _minpoly_compose(ex, x, dom)
    if not pw.is_rational:
        raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)
    if pw < 0:
        if mp == x:
            raise ZeroDivisionError('%s is zero' % ex)
        mp = _invertx(mp, x)
        if pw == -1:
            return mp
        pw = -pw
        ex = 1/ex

    y = Dummy(str(x))
    mp = mp.subs({x: y})
    n, d = pw.as_numer_denom()
    res = Poly(resultant(mp, x**d - y**n, gens=[y]), x, domain=dom)
    _, factors = res.factor_list()
    res = _choose_factor(factors, x, ex**pw, dom)
    return res.as_expr()

    def __new__(cls, expr, coeffs=Tuple(), alias=None, **args):
        """Construct a new algebraic number. """
        expr = sympify(expr)

        if isinstance(expr, (tuple, Tuple)):
            minpoly, root = expr

            if not minpoly.is_Poly:
                minpoly = Poly(minpoly)
        elif expr.is_AlgebraicNumber:
            minpoly, root = expr.minpoly, expr.root
        else:
            minpoly, root = minimal_polynomial(
                expr, args.get('gen'), polys=True), expr

        dom = minpoly.get_domain()

        if coeffs != Tuple():
            if not isinstance(coeffs, ANP):
                rep = DMP.from_sympy_list(sympify(coeffs), 0, dom)
                scoeffs = Tuple(*coeffs)
            else:
                rep = DMP.from_list(coeffs.to_list(), 0, dom)
                scoeffs = Tuple(*coeffs.to_list())

            if rep.degree() >= minpoly.degree():
                rep = rep.rem(minpoly.rep)

            sargs = (root, scoeffs)

        else:
            rep = DMP.from_list([1, 0], 0, dom)

            if ask(Q.negative(root)):
                rep = -rep

            sargs = (root, coeffs)

        if alias is not None:
            if not isinstance(alias, Symbol):
                alias = Symbol(alias)
            sargs = sargs + (alias,)

        obj = Expr.__new__(cls, *sargs)

        obj.rep = rep
        obj.root = root
        obj.alias = alias
        obj.minpoly = minpoly

        return obj

def full_coeffs(p,D,s):
    '''
    Computes coefficients of a polynomial matrix
    p : polynomial
    D :degree  up to which we want to complete the list 	with zeroes
    s: variable of the polynomial matrix
    Example:  TODO
    '''
    p=Poly(p,s,domain='QQ')  # in order to use the all_coeffs method of polynomial class
    c=p.all_coeffs()
    if len(c)==D+1:
        return c
    else:
        difference=D+1-len(c)
        return difference*[0] +c

def chebyshevt_poly(n, x=None, **args):
    """Generates Chebyshev polynomial of the first kind of degree `n` in `x`. """
    if n < 0:
        raise ValueError("can't generate 1st kind Chebyshev polynomial of degree %s" % n)

    if x is not None:
        x = sympify(x)
    else:
        x = Symbol('x', dummy=True)

    poly = Poly(DMP(dup_chebyshevt(int(n), ZZ), ZZ), x)

    if not args.get('polys', False):
        return poly.as_basic()
    else:
        return poly

def cyclotomic_poly(n, x=None, **args):
    """Generates cyclotomic polynomial of order `n` in `x`. """
    if n <= 0:
        raise ValueError("can't generate cyclotomic polynomial of order %s" % n)

    if x is not None:
        x = sympify(x)
    else:
        x = Symbol('x', dummy=True)

    poly = Poly(DMP(dup_zz_cyclotomic_poly(int(n), ZZ), ZZ), x)

    if not args.get('polys', False):
        return poly.as_basic()
    else:
        return poly

def gegenbauer_poly(n, a, x=None, polys=False):
    """Generates Gegenbauer polynomial of degree `n` in `x`.

    Parameters
    ==========

    n : int
        `n` decides the degree of polynomial
    x : optional
    a
        Decides minimal domain for the list of
        coefficients.
    polys : bool, optional
        ``polys=True`` returns an expression, otherwise
        (default) returns an expression.
    """
    if n < 0:
        raise ValueError(
            "can't generate Gegenbauer polynomial of degree %s" % n)

    K, a = construct_domain(a, field=True)
    poly = DMP(dup_gegenbauer(int(n), a, K), K)

    if x is not None:
        poly = Poly.new(poly, x)
    else:
        poly = PurePoly.new(poly, Dummy('x'))

    return poly if polys else poly.as_expr()

def chebyshevu_poly(n, x=None, **args):
    """Generates Chebyshev polynomial of the second kind of degree `n` in `x`. """
    if n < 0:
        raise ValueError("can't generate 2nd kind Chebyshev polynomial of degree %s" % n)

    if x is not None:
        x = sympify(x)
    else:
        x = Dummy('x')

    poly = Poly(DMP(dup_chebyshevu(int(n), ZZ), ZZ), x)

    if not args.get('polys', False):
        return poly.as_basic()
    else:
        return poly

def jacobi_poly(n, a, b, x=None, polys=False):
    """Generates Jacobi polynomial of degree `n` in `x`.

    Parameters
    ==========

    n : int
        `n` decides the degree of polynomial
    a
        Lower limit of minimal domain for the list of
        coefficients.
    b
        Upper limit of minimal domain for the list of
        coefficients.
    x : optional
    polys : bool, optional
        ``polys=True`` returns an expression, otherwise
        (default) returns an expression.
    """
    if n < 0:
        raise ValueError("can't generate Jacobi polynomial of degree %s" % n)

    K, v = construct_domain([a, b], field=True)
    poly = DMP(dup_jacobi(int(n), v[0], v[1], K), K)

    if x is not None:
        poly = Poly.new(poly, x)
    else:
        poly = PurePoly.new(poly, Dummy('x'))

    return poly if polys else poly.as_expr()

def laguerre_poly(n, x=None, alpha=None, polys=False):
    """Generates Laguerre polynomial of degree `n` in `x`.

    Parameters
    ==========

    n : int
        `n` decides the degree of polynomial
    x : optional
    alpha
        Decides minimal domain for the list
        of coefficients.
    polys : bool, optional
        ``polys=True`` returns an expression, otherwise
        (default) returns an expression.
    """
    if n < 0:
        raise ValueError("can't generate Laguerre polynomial of degree %s" % n)

    if alpha is not None:
        K, alpha = construct_domain(
            alpha, field=True)  # XXX: ground_field=True
    else:
        K, alpha = QQ, QQ(0)

    poly = DMP(dup_laguerre(int(n), alpha, K), K)

    if x is not None:
        poly = Poly.new(poly, x)
    else:
        poly = PurePoly.new(poly, Dummy('x'))

    return poly if polys else poly.as_expr()

def chebyshevu_poly(n, x=None, polys=False):
    """Generates Chebyshev polynomial of the second kind of degree `n` in `x`.

    Parameters
    ==========

    n : int
        `n` decides the degree of polynomial
    x : optional
    polys : bool, optional
        ``polys=True`` returns an expression, otherwise
        (default) returns an expression.
    """
    if n < 0:
        raise ValueError(
            "can't generate 2nd kind Chebyshev polynomial of degree %s" % n)

    poly = DMP(dup_chebyshevu(int(n), ZZ), ZZ)

    if x is not None:
        poly = Poly.new(poly, x)
    else:
        poly = PurePoly.new(poly, Dummy('x'))

    return poly if polys else poly.as_expr()

def hermite_poly(n, x=None, **args):
    """Generates Hermite polynomial of degree `n` in `x`. """
    if n < 0:
        raise ValueError("can't generate Hermite polynomial of degree %s" % n)

    if x is not None:
        x = sympify(x)
    else:
        x = Dummy('x')

    poly = Poly(DMP(dup_hermite(int(n), ZZ), ZZ), x)

    if not args.get('polys', False):
        return poly.as_basic()
    else:
        return poly

def laguerre_poly(n, x=None, **args):
    """Generates Laguerre polynomial of degree `n` in `x`. """
    if n < 0:
        raise ValueError("can't generate Laguerre polynomial of degree %s" % n)

    if x is not None:
        x = sympify(x)
    else:
        x = Symbol('x', dummy=True)

    poly = Poly(DMP(dup_laguerre(int(n), QQ), QQ), x)

    if not args.get('polys', False):
        return poly.as_basic()
    else:
        return poly

    def _try_decompose(f):
        """Find roots using functional decomposition. """
        
        factors, roots = f.decompose(), []
        if len(factors) == 1:
            for root in _try_heuristics(factors[0]):
                roots.append(root)
            return (roots, False)

        add_comment("Use the substitution")
        t = Dummy("t1")
        add_eq(t, compose_(factors[1:]).as_expr())
        g = Poly(factors[0].as_expr().subs(f.gen, t), t)
        add_comment("We have")
        add_eq(g.as_expr(), 0)
        for root in _try_heuristics(g):
            roots.append(root)

        # f(x) = f1(f2(f3(x))) = 0
        # f(x) = f1(t) = 0 --> f2(f3(x)) = t1, t2, ..., tn

        i = 1
        for factor in factors[1:]:
            previous, roots = list(roots), []

            h = compose_(factors[i:]).as_expr()
            add_comment("Therefore")
            for root in previous:
                add_eq(h, root)

            for root in previous:
                if i < len(factors) - 1:
                    add_comment("Solve the equation")
                    add_eq(h, root)
                    add_comment("Use the substitution")
                    t = Dummy("t" + str(i+1))
                    add_eq(t, compose_(factors[(i+1):]).as_expr())
                    g = Poly(factor.as_expr().subs(f.gen, t), t) - Poly(root, t)
                else:
                    g = factor - Poly(root, f.gen)

                for root in _try_heuristics(g):
                    roots.append(root)
            i += 1

        return (roots, True)