def test_positive():
    x, y, z, w = symbols('xyzw')
    assert ask(x, Q.positive, Assume(x, Q.positive)) == True
    assert ask(x, Q.positive, Assume(x, Q.negative)) == False
    assert ask(x, Q.positive, Assume(x, Q.nonzero)) == None

    assert ask(-x, Q.positive, Assume(x, Q.positive)) == False
    assert ask(-x, Q.positive, Assume(x, Q.negative)) == True

    assert ask(x+y, Q.positive, Assume(x, Q.positive) & \
                     Assume(y, Q.positive)) == True
    assert ask(x+y, Q.positive, Assume(x, Q.positive) & \
                     Assume(y, Q.negative)) == None

    assert ask(2*x,  Q.positive, Assume(x, Q.positive)) == True
    assumptions =  Assume(x, Q.positive) & Assume(y, Q.negative) & \
                    Assume(z, Q.negative) & Assume(w, Q.positive)
    assert ask(x*y*z,  Q.positive)  == None
    assert ask(x*y*z,  Q.positive, assumptions) == True
    assert ask(-x*y*z, Q.positive, assumptions) == False

    assert ask(x**2, Q.positive, Assume(x, Q.positive)) == True
    assert ask(x**2, Q.positive, Assume(x, Q.negative)) == True

    #exponential
    assert ask(exp(x),     Q.positive, Assume(x, Q.real)) == True
    assert ask(x + exp(x), Q.positive, Assume(x, Q.real)) == None

    #absolute value
    assert ask(Abs(x), Q.positive) == None # Abs(0) = 0
    assert ask(Abs(x), Q.positive, Assume(x, Q.positive)) == True

def moveup2(s, x):
    r = SubsSet()
    for expr, var in s.iteritems():
        r[expr.subs(x, exp(x))] = var
    for var, expr in s.rewrites.iteritems():
        r.rewrites[var] = s.rewrites[var].subs(x, exp(x))
    return r

def test_contraction_structure_Mul_and_Pow():
    x = IndexedBase('x')
    y = IndexedBase('y')
    i, j, k = Idx('i'), Idx('j'), Idx('k')

    i_ji = x[i]**(y[j]*x[i])
    assert get_contraction_structure(i_ji) == {None: set([i_ji])}
    ij_i = (x[i]*y[j])**(y[i])
    assert get_contraction_structure(ij_i) == {None: set([ij_i])}
    j_ij_i = x[j]*(x[i]*y[j])**(y[i])
    assert get_contraction_structure(j_ij_i) == {(j,): set([j_ij_i])}
    j_i_ji = x[j]*x[i]**(y[j]*x[i])
    assert get_contraction_structure(j_i_ji) == {(j,): set([j_i_ji])}
    ij_exp_kki = x[i]*y[j]*exp(y[i]*y[k, k])
    result = get_contraction_structure(ij_exp_kki)
    expected = {
        (i,): set([ij_exp_kki]),
        ij_exp_kki: [{
                     None: set([exp(y[i]*y[k, k])]),
                exp(y[i]*y[k, k]): [{
                    None: set([y[i]*y[k, k]]),
                    y[i]*y[k, k]: [{(k,): set([y[k, k]])}]
                }]}
        ]
    }
    assert result == expected

def test_positive():
    x, y, z, w = symbols('x,y,z,w')
    assert ask(Q.positive(x), Q.positive(x)) == True
    assert ask(Q.positive(x), Q.negative(x)) == False
    assert ask(Q.positive(x), Q.nonzero(x)) == None

    assert ask(Q.positive(-x), Q.positive(x)) == False
    assert ask(Q.positive(-x), Q.negative(x)) == True

    assert ask(Q.positive(x+y), Q.positive(x) & Q.positive(y)) == True
    assert ask(Q.positive(x+y), Q.positive(x) & Q.negative(y)) == None

    assert ask(Q.positive(2*x), Q.positive(x)) == True
    assumptions =  Q.positive(x) & Q.negative(y) & Q.negative(z) & Q.positive(w)
    assert ask(Q.positive(x*y*z))  == None
    assert ask(Q.positive(x*y*z), assumptions) == True
    assert ask(Q.positive(-x*y*z), assumptions) == False

    assert ask(Q.positive(x**2), Q.positive(x)) == True
    assert ask(Q.positive(x**2), Q.negative(x)) == True

    #exponential
    assert ask(Q.positive(exp(x)), Q.real(x)) == True
    assert ask(Q.positive(x + exp(x)), Q.real(x)) == None

    #absolute value
    assert ask(Q.positive(Abs(x))) == None # Abs(0) = 0
    assert ask(Q.positive(Abs(x)), Q.positive(x)) == True

def _generate_patterns():
    """
    Generates patterns for transcendental equations.

    This is lazily calculated (called) in the tsolve() function and stored in
    the patterns global variable.
    """

    tmp1 = _f ** (_h-(_c*_g/_b))
    tmp2 = (-_e*tmp1/_a)**(1/_d)
    global _patterns
    _patterns = [
        (_a*(_b*_x+_c)**_d + _e   ,
            ((-(_e/_a))**(1/_d)-_c)/_b),
        (_b+_c*exp(_d*_x+_e) ,
            (log(-_b/_c)-_e)/_d),
        (_a*_x+_b+_c*exp(_d*_x+_e) ,
            -_b/_a-LambertW(_c*_d*exp(_e-_b*_d/_a)/_a)/_d),
        (_b+_c*_f**(_d*_x+_e) ,
            (log(-_b/_c)-_e*log(_f))/_d/log(_f)),
        (_a*_x+_b+_c*_f**(_d*_x+_e) ,
            -_b/_a-LambertW(_c*_d*_f**(_e-_b*_d/_a)*log(_f)/_a)/_d/log(_f)),
        (_b+_c*log(_d*_x+_e) ,
            (exp(-_b/_c)-_e)/_d),
        (_a*_x+_b+_c*log(_d*_x+_e) ,
            -_e/_d+_c/_a*LambertW(_a/_c/_d*exp(-_b/_c+_a*_e/_c/_d))),
        (_a*(_b*_x+_c)**_d + _e*_f**(_g*_x+_h) ,
            -_c/_b-_d*LambertW(-tmp2*_g*log(_f)/_b/_d)/_g/log(_f))
    ]

 def _expr_big_minus(cls, a, z, n):
     from sympy import I, pi, exp, sqrt, atan, sin
     if n.is_even:
         return (1 + z)**a*exp(2*pi*I*n*a)*sqrt(z)*sin(2*a*atan(sqrt(z)))
     else:
         return (1 + z)**a*exp(2*pi*I*n*a)*sqrt(z) \
                *sin(2*a*atan(sqrt(z)) - 2*pi*a)

def mrv(e, x):
    """Returns a SubsSet of most rapidly varying (mrv) subexpressions of 'e',
       and e rewritten in terms of these"""
    e = powsimp(e, deep=True, combine='exp')
    assert isinstance(e, Basic)
    if not e.has(x):
        return SubsSet(), e
    elif e == x:
        s = SubsSet()
        return s, s[x]
    elif e.is_Mul or e.is_Add:
        i, d = e.as_independent(x)  # throw away x-independent terms
        if d.func != e.func:
            s, expr = mrv(d, x)
            return s, e.func(i, expr)
        a, b = d.as_two_terms()
        s1, e1 = mrv(a, x)
        s2, e2 = mrv(b, x)
        return mrv_max1(s1, s2, e.func(i, e1, e2), x)
    elif e.is_Pow:
        b, e = e.as_base_exp()
        if e.has(x):
            return mrv(exp(e * log(b)), x)
        else:
            s, expr = mrv(b, x)
            return s, expr**e
    elif e.func is log:
        s, expr = mrv(e.args[0], x)
        return s, log(expr)
    elif e.func is exp:
        # We know from the theory of this algorithm that exp(log(...)) may always
        # be simplified here, and doing so is vital for termination.
        if e.args[0].func is log:
            return mrv(e.args[0].args[0], x)
        if limitinf(e.args[0], x).is_unbounded:
            s1 = SubsSet()
            e1 = s1[e]
            s2, e2 = mrv(e.args[0], x)
            su = s1.union(s2)[0]
            su.rewrites[e1] = exp(e2)
            return mrv_max3(s1, e1, s2, exp(e2), su, e1, x)
        else:
            s, expr = mrv(e.args[0], x)
            return s, exp(expr)
    elif e.is_Function:
        l = [mrv(a, x) for a in e.args]
        l2 = [s for (s, _) in l if s != SubsSet()]
        if len(l2) != 1:
            # e.g. something like BesselJ(x, x)
            raise NotImplementedError("MRV set computation for functions in"
                                      " several variables not implemented.")
        s, ss = l2[0], SubsSet()
        args = [ss.do_subs(x[1]) for x in l]
        return s, e.func(*args)
    elif e.is_Derivative:
        raise NotImplementedError("MRV set computation for derviatives"
                                  " not implemented yet.")
        return mrv(e.args[0], x)
    raise NotImplementedError(
        "Don't know how to calculate the mrv of '%s'" % e)

def moveup2(s, x):
    r = SubsSet()
    for expr, var in s.items():
        r[expr.xreplace({x: exp(x)})] = var
    for var, expr in s.rewrites.items():
        r.rewrites[var] = s.rewrites[var].xreplace({x: exp(x)})
    return r

 def _expr_big(cls, a, z, n):
     if n.is_even:
         return ((sqrt(z) + 1)**(2*a)*exp(2*pi*I*n*a) +
                 (sqrt(z) - 1)**(2*a)*exp(2*pi*I*(n - 1)*a))/2
     else:
         n -= 1
         return ((sqrt(z) - 1)**(2*a)*exp(2*pi*I*a*(n + 1)) +
                 (sqrt(z) + 1)**(2*a)*exp(2*pi*I*a*n))/2

def test_wavefunction():
  Psi = {
    0: (nu/pi)**(S(1)/4) * exp(-nu * x**2 /2),
    1: (nu/pi)**(S(1)/4) * sqrt(2*nu) * x * exp(-nu * x**2 /2),
    2: (nu/pi)**(S(1)/4) * (2 * nu * x**2 - 1)/sqrt(2) * exp(-nu * x**2 /2),
    3: (nu/pi)**(S(1)/4) * sqrt(nu/3) * (2 * nu * x**3 - 3 * x) * exp(-nu * x**2 /2)
  }
  for n in Psi:
    assert simplify(psi_n(n, x, m, omega) - Psi[n]) == 0

 def _expr_big(cls, a, z, n):
     from sympy import I, pi, exp, sqrt, atan, sin
     if n.is_even:
         return sqrt(z)/2*((sqrt(z) - 1)**(2*a)*exp(2*pi*I*a*(n - 1)) -
                           (sqrt(z) + 1)**(2*a)*exp(2*pi*I*a*n))
     else:
         n -= 1
         return sqrt(z)/2*((sqrt(z) - 1)**(2*a)*exp(2*pi*I*a*(n + 1)) -
                           (sqrt(z) + 1)**(2*a)*exp(2*pi*I*a*n))

def test_exp():
    A = MatrixSymbol('A', 2, 2)
    B = MatrixSymbol('B', 2, 2)
    expr1 = exp(A)*exp(B)
    expr2 = exp(B)*exp(A)
    assert expr1 != expr2
    assert expr1 - expr2 != 0
    assert not isinstance(expr1, exp)
    assert not isinstance(expr2, exp)

def rs_exp(p, x, prec):
    """
    Exponentiation of a series modulo ``O(x**prec)``

    Examples
    ========

    >>> from sympy.polys.domains import QQ
    >>> from sympy.polys.rings import ring
    >>> from sympy.polys.ring_series import rs_exp
    >>> R, x = ring('x', QQ)
    >>> rs_exp(x**2, x, 7)
    1/6*x**6 + 1/2*x**4 + x**2 + 1
    """
    if rs_is_puiseux(p, x):
        return rs_puiseux(rs_exp, p, x, prec)
    R = p.ring
    c = _get_constant_term(p, x)
    if c:
        if R.domain is EX:
            c_expr = c.as_expr()
            const = exp(c_expr)
        elif isinstance(c, PolyElement):
            try:
                c_expr = c.as_expr()
                const = R(exp(c_expr))
            except ValueError:
                R = R.add_gens([exp(c_expr)])
                p = p.set_ring(R)
                x = x.set_ring(R)
                c = c.set_ring(R)
                const = R(exp(c_expr))
        else:
            try:
                const = R(exp(c))
            except ValueError:
                    raise DomainError("The given series can't be expanded in "
                                      "this domain.")
        p1 = p - c

    # Makes use of sympy fuctions to evaluate the values of the cos/sin
    # of the constant term.
        return const*rs_exp(p1, x, prec)

    if len(p) > 20:
        return _exp1(p, x, prec)
    one = R(1)
    n = 1
    k = 1
    c = []
    for k in range(prec):
        c.append(one/n)
        k += 1
        n *= k

    r = rs_series_from_list(p, c, x, prec)
    return r

def solve_ODE_second_order(eq, f):
    """
    solves many kinds of second order odes, different methods are used
    depending on the form of the given equation. So far the constants
    coefficients case and a special case are implemented.
    """
    x = f.args[0]
    f = f.func

    #constant coefficients case: af''(x)+bf'(x)+cf(x)=0
    a = Wild('a', exclude=[x])
    b = Wild('b', exclude=[x])
    c = Wild('c', exclude=[x])

    r = eq.match(a*f(x).diff(x,x) + c*f(x))
    if r:
        return Symbol("C1")*C.sin(sqrt(r[c]/r[a])*x)+Symbol("C2")*C.cos(sqrt(r[c]/r[a])*x)

    r = eq.match(a*f(x).diff(x,x) + b*diff(f(x),x) + c*f(x))
    if r:
        r1 = solve(r[a]*x**2 + r[b]*x + r[c], x)
        if r1[0].is_real:
            if len(r1) == 1:
                return (Symbol("C1") + Symbol("C2")*x)*exp(r1[0]*x)
            else:
                return Symbol("C1")*exp(r1[0]*x) + Symbol("C2")*exp(r1[1]*x)
        else:
            r2 = abs((r1[0] - r1[1])/(2*S.ImaginaryUnit))
            return (Symbol("C2")*C.cos(r2*x) + Symbol("C1")*C.sin(r2*x))*exp((r1[0] + r1[1])*x/2)

    #other cases of the second order odes will be implemented here

    #special equations, that we know how to solve
    a = Wild('a')
    t = x*C.exp(f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        return -solve_ODE_1(f(x), x)

    t = x*C.exp(-f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( r[a]*t.diff(x,2)/t ) == eq.subs(f, t)
        return solve_ODE_1(f(x), x)

    neq = eq*C.exp(f(x))/C.exp(-f(x))
    r = neq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( t.diff(x,2)*r[a]/t ).expand() == eq
        return solve_ODE_1(f(x), x)

    raise NotImplementedError("solve_ODE_second_order: cannot solve " + str(eq))

def test_functional_exponent():
    t = standard_transformations + (convert_xor, function_exponentiation)
    x = Symbol('x')
    y = Symbol('y')
    a = Symbol('a')
    yfcn = Function('y')
    assert parse_expr("sin^2(x)", transformations=t) == (sin(x))**2
    assert parse_expr("sin^y(x)", transformations=t) == (sin(x))**y
    assert parse_expr("exp^y(x)", transformations=t) == (exp(x))**y
    assert parse_expr("E^y(x)", transformations=t) == exp(yfcn(x))
    assert parse_expr("a^y(x)", transformations=t) == a**(yfcn(x))

def test_rs_series():
    x, a, b, c = symbols('x, a, b, c')

    assert rs_series(a, a, 5).as_expr() == a
    assert rs_series(sin(a), a, 5).as_expr() == (sin(a).series(a, 0,
        5)).removeO()
    assert rs_series(sin(a) + cos(a), a, 5).as_expr() == ((sin(a) +
        cos(a)).series(a, 0, 5)).removeO()
    assert rs_series(sin(a)*cos(a), a, 5).as_expr() == ((sin(a)*
        cos(a)).series(a, 0, 5)).removeO()

    p = (sin(a) - a)*(cos(a**2) + a**4/2)
    assert expand(rs_series(p, a, 10).as_expr()) == expand(p.series(a, 0,
        10).removeO())

    p = sin(a**2/2 + a/3) + cos(a/5)*sin(a/2)**3
    assert expand(rs_series(p, a, 5).as_expr()) == expand(p.series(a, 0,
        5).removeO())

    p = sin(x**2 + a)*(cos(x**3 - 1) - a - a**2)
    assert expand(rs_series(p, a, 5).as_expr()) == expand(p.series(a, 0,
        5).removeO())

    p = sin(a**2 - a/3 + 2)**5*exp(a**3 - a/2)
    assert expand(rs_series(p, a, 10).as_expr()) == expand(p.series(a, 0,
        10).removeO())

    p = sin(a + b + c)
    assert expand(rs_series(p, a, 5).as_expr()) == expand(p.series(a, 0,
        5).removeO())

    p = tan(sin(a**2 + 4) + b + c)
    assert expand(rs_series(p, a, 6).as_expr()) == expand(p.series(a, 0,
        6).removeO())

    p = a**QQ(2,5) + a**QQ(2,3) + a

    r = rs_series(tan(p), a, 2)
    assert r.as_expr() == a**QQ(9,5) + a**QQ(26,15) + a**QQ(22,15) + a**QQ(6,5)/3 + \
        a + a**QQ(2,3) + a**QQ(2,5)

    r = rs_series(exp(p), a, 1)
    assert r.as_expr() == a**QQ(4,5)/2 + a**QQ(2,3) + a**QQ(2,5) + 1

    r = rs_series(sin(p), a, 2)
    assert r.as_expr() == -a**QQ(9,5)/2 - a**QQ(26,15)/2 - a**QQ(22,15)/2 - \
        a**QQ(6,5)/6 + a + a**QQ(2,3) + a**QQ(2,5)

    r = rs_series(cos(p), a, 2)
    assert r.as_expr() == a**QQ(28,15)/6 - a**QQ(5,3) + a**QQ(8,5)/24 - a**QQ(7,5) - \
        a**QQ(4,3)/2 - a**QQ(16,15) - a**QQ(4,5)/2 + 1

    assert rs_series(sin(a)/7, a, 5).as_expr() == (sin(a)/7).series(a, 0,
            5).removeO()

def test_bounded():
    x, y = symbols('x,y')
    assert ask(Q.bounded(x)) == False
    assert ask(Q.bounded(x), Q.bounded(x)) == True
    assert ask(Q.bounded(x), Q.bounded(y)) == False
    assert ask(Q.bounded(x), Q.complex(x)) == False

    assert ask(Q.bounded(x+1)) == False
    assert ask(Q.bounded(x+1), Q.bounded(x)) == True
    assert ask(Q.bounded(x+y)) == None
    assert ask(Q.bounded(x+y), Q.bounded(x)) == False
    assert ask(Q.bounded(x+1), Q.bounded(x) & Q.bounded(y)) == True

    assert ask(Q.bounded(2*x)) == False
    assert ask(Q.bounded(2*x), Q.bounded(x)) == True
    assert ask(Q.bounded(x*y)) == None
    assert ask(Q.bounded(x*y), Q.bounded(x)) == False
    assert ask(Q.bounded(x*y), Q.bounded(x) & Q.bounded(y)) == True

    assert ask(Q.bounded(x**2)) == False
    assert ask(Q.bounded(2**x)) == False
    assert ask(Q.bounded(2**x), Q.bounded(x)) == True
    assert ask(Q.bounded(x**x)) == False
    assert ask(Q.bounded(Rational(1,2) ** x)) == None
    assert ask(Q.bounded(Rational(1,2) ** x), Q.positive(x)) == True
    assert ask(Q.bounded(Rational(1,2) ** x), Q.negative(x)) == False
    assert ask(Q.bounded(sqrt(x))) == False

    # sign function
    assert ask(Q.bounded(sign(x))) == True
    assert ask(Q.bounded(sign(x)), ~Q.bounded(x)) == True

    # exponential functions
    assert ask(Q.bounded(log(x))) == False
    assert ask(Q.bounded(log(x)), Q.bounded(x)) == True
    assert ask(Q.bounded(exp(x))) == False
    assert ask(Q.bounded(exp(x)), Q.bounded(x)) == True
    assert ask(Q.bounded(exp(2))) == True

    # trigonometric functions
    assert ask(Q.bounded(sin(x))) == True
    assert ask(Q.bounded(sin(x)), ~Q.bounded(x)) == True
    assert ask(Q.bounded(cos(x))) == True
    assert ask(Q.bounded(cos(x)), ~Q.bounded(x)) == True
    assert ask(Q.bounded(2*sin(x))) == True
    assert ask(Q.bounded(sin(x)**2)) == True
    assert ask(Q.bounded(cos(x)**2)) == True
    assert ask(Q.bounded(cos(x) + sin(x))) == True

def test_bounded():
    x, y = symbols('xy')
    assert ask(x, Q.bounded) == False
    assert ask(x, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x, Q.bounded, Assume(y, Q.bounded)) == False
    assert ask(x, Q.bounded, Assume(x, Q.complex)) == False

    assert ask(x+1, Q.bounded) == False
    assert ask(x+1, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x+y, Q.bounded) == None
    assert ask(x+y, Q.bounded, Assume(x, Q.bounded)) == False
    assert ask(x+1, Q.bounded, Assume(x, Q.bounded) & \
                Assume(y, Q.bounded)) == True

    assert ask(2*x, Q.bounded) == False
    assert ask(2*x, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x*y, Q.bounded) == None
    assert ask(x*y, Q.bounded, Assume(x, Q.bounded)) == False
    assert ask(x*y, Q.bounded, Assume(x, Q.bounded) & \
                Assume(y, Q.bounded)) == True

    assert ask(x**2, Q.bounded) == False
    assert ask(2**x, Q.bounded) == False
    assert ask(2**x, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x**x, Q.bounded) == False
    assert ask(Rational(1,2) ** x, Q.bounded) == True
    assert ask(x ** Rational(1,2), Q.bounded) == False

    # sign function
    assert ask(sign(x), Q.bounded) == True
    assert ask(sign(x), Q.bounded, Assume(x, Q.bounded, False)) == True

    # exponential functions
    assert ask(log(x), Q.bounded) == False
    assert ask(log(x), Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(exp(x), Q.bounded) == False
    assert ask(exp(x), Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(exp(2), Q.bounded) == True

    # trigonometric functions
    assert ask(sin(x), Q.bounded) == True
    assert ask(sin(x), Q.bounded, Assume(x, Q.bounded, False)) == True
    assert ask(cos(x), Q.bounded) == True
    assert ask(cos(x), Q.bounded, Assume(x, Q.bounded, False)) == True
    assert ask(2*sin(x), Q.bounded) == True
    assert ask(sin(x)**2, Q.bounded) == True
    assert ask(cos(x)**2, Q.bounded) == True
    assert ask(cos(x) + sin(x), Q.bounded) == True

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 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)