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 is_primitive_root(a, p):
    """
    Returns True if ``a`` is a primitive root of ``p``

    ``a`` is said to be the primitive root of ``p`` if gcd(a, p) == 1 and
    totient(p) is the smallest positive number s.t.

        a**totient(p) cong 1 mod(p)

    Examples
    ========

    >>> from sympy.ntheory import is_primitive_root, n_order, totient
    >>> is_primitive_root(3, 10)
    True
    >>> is_primitive_root(9, 10)
    False
    >>> n_order(3, 10) == totient(10)
    True
    >>> n_order(9, 10) == totient(10)
    False

    """
    a, p = as_int(a), as_int(p)
    if igcd(a, p) != 1:
        raise ValueError("The two numbers should be relatively prime")
    if a > p:
        a = a % p
    return n_order(a, p) == totient(p)

def is_primitive_root(a, p):
    """
    Returns True if ``a`` is a primitive root of ``n``

    ``a`` is said to be the primitive root of ``n`` if gcd(a, n) == 1 and
    totient(n) is the smallest positive number s.t.

        a**totient(n) cong 1 mod(n)

    **Examples**
    >>> from sympy.ntheory import is_primitive_root, n_order, totient
    >>> is_primitive_root(3, 10)
    True
    >>> is_primitive_root(9, 10)
    False
    >>> n_order(3, 10) == totient(10)
    True
    >>> n_order(9, 10) == totient(10)
    False
    """
    a, p = int_tested(a, p)
    if igcd(a, p) != 1:
        raise ValueError("The two numbers should be relatively prime")
    if a > p:
        a = a % p
    if n_order(a, p) == totient_(p):
        return True
    else:
        return False

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 deflate(f, *G):
        ring = f.ring
        polys = [f] + list(G)

        J = [0]*ring.ngens

        for p in polys:
            for monom in p.monoms():
                for i, m in enumerate(monom):
                    J[i] = igcd(J[i], m)

        for i, b in enumerate(J):
            if not b:
                J[i] = 1

        J = tuple(J)

        if all(b == 1 for b in J):
            return J, polys

        H = []

        for p in polys:
            h = ring.zero

            for I, coeff in p.terms():
                N = [ i // j for i, j in zip(I, J) ]
                h[tuple(N)] = coeff

            H.append(h)

        return J, H

def n_order(a, n):
    """Returns the order of ``a`` modulo ``n``.

    The order of ``a`` modulo ``n`` is the smallest integer
    ``k`` such that ``a**k`` leaves a remainder of 1 with ``n``.

    Examples
    ========

    >>> from sympy.ntheory import n_order
    >>> n_order(3, 7)
    6
    >>> n_order(4, 7)
    3
    """
    a, n = int_tested(a, n)
    if igcd(a, n) != 1:
        raise ValueError("The two numbers should be relatively prime")
    group_order = totient(n)
    factors = factorint(group_order)
    order = 1
    if a > n:
        a = a % n
    for p, e in factors.iteritems():
        exponent = group_order
        for f in xrange(0, e + 1):
            if (a ** (exponent)) % n != 1:
                order *= p ** (e - f + 1)
                break
            exponent = exponent // p
    return order

def poly_reduce(f, g, *symbols):
    """Removes common content from a pair of polynomials.

       >>> from sympy import *
       >>> x = Symbol('x')

       >>> f = Poly(2930944*x**6 + 2198208*x**4 + 549552*x**2 + 45796, x)
       >>> g = Poly(17585664*x**5 + 8792832*x**3 + 1099104*x, x)

       >>> F, G = poly_reduce(f, g)

       >>> F
       Poly(64*x**6 + 48*x**4 + 12*x**2 + 1, x)
       >>> G
       Poly(384*x**5 + 192*x**3 + 24*x, x)

    """
    if not isinstance(f, Poly):
        f = Poly(f, *symbols)
    elif symbols:
        raise SymbolsError("Redundant symbols were given")

    f, g = f.unify_with(g)

    fc = int(f.content)
    gc = int(g.content)

    cont = igcd(fc, gc)

    if cont != 1:
        f = f.div_term(cont)
        g = g.div_term(cont)

    return f, g

def _a(n, j, prec):
    """Compute the inner sum in the HRR formula."""
    if j == 1:
        return fone
    s = fzero
    pi = pi_fixed(prec)
    for h in xrange(1, j):
        if igcd(h, j) != 1:
            continue
        # & with mask to compute fractional part of fixed-point number
        one = 1 << prec
        onemask = one - 1
        half = one >> 1
        g = 0
        if j >= 3:
            for k in xrange(1, j):
                t = h*k*one//j
                if t > 0:
                    frac = t & onemask
                else:
                    frac = -((-t) & onemask)
                g += k*(frac - half)
        g = ((g - 2*h*n*one)*pi//j) >> prec
        s = mpf_add(s, mpf_cos(from_man_exp(g, -prec), prec), prec)
    return s

def pollard_pm1(n, B=10, seed=1234):
    """Use Pollard's p-1 method to try to extract a factor of n. The
    returned factor may be a composite number. The search is performed
    up to a smoothness bound B; if no factor is found, None is
    returned.

    The p-1 algorithm is a Monte Carlo method whose outcome can
    be affected by changing the random seed value.

    Example usage
    =============
    With the default smoothness bound, this number can't be cracked:
        >>> pollard_pm1(21477639576571)

    Increasing the smoothness bound helps:
        >>> pollard_pm1(21477639576571, 2000)
        4410317

    References
    ==========
    Richard Crandall & Carl Pomerance (2005), "Prime Numbers:
    A Computational Perspective", Springer, 2nd edition, 236-238

    """
    from math import log
    prng = random.Random(seed + B)
    a = prng.randint(2, n-1)
    for p in sieve.primerange(2, B):
        e = int(log(B, p))
        a = pow(a, p**e, n)
    g = igcd(a-1, n)
    if 1 < g < n:
        return g
    else:
        return None

def pollard_rho(n, max_iters=5, seed=1234):
    """Use Pollard's rho method to try to extract a factor of n. The
    returned factor may be a composite number. A maximum of max_iters
    iterations are performed; if no factor is found, None is returned.

    The rho algorithm is a Monte Carlo method whose outcome can
    be affected by changing the random seed value.

    References
    ==========
    Richard Crandall & Carl Pomerance (2005), "Prime Numbers:
    A Computational Perspective", Springer, 2nd edition, 229-231

    """
    prng = random.Random(seed + max_iters)
    for i in range(max_iters):
        # Alternative good nonrandom choice: a = 1
        a = prng.randint(1, n-3)
        # Alternative good nonrandom choice: s = 2
        s = prng.randint(0, n-1)
        U = V = s
        F = lambda x: (x**2 + a) % n
        while 1:
            U = F(U)
            V = F(F(V))
            g = igcd(abs(U-V), n)
            if g == 1:
                continue
            if g == n:
                break
            return g
    return None

def jacobi_symbol(m, n):
    """
    Returns the product of the legendre_symbol(m, p)
    for all the prime factors p of n.

    Returns
    =======

    1. 0 if m cong 0 mod(n)
    2. 1 if x**2 cong m mod(n) has a solution
    3. -1 otherwise

    Examples
    ========

    >>> from sympy.ntheory import jacobi_symbol, legendre_symbol
    >>> from sympy import Mul, S
    >>> jacobi_symbol(45, 77)
    -1
    >>> jacobi_symbol(60, 121)
    1

    The relationship between the jacobi_symbol and legendre_symbol can
    be demonstrated as follows:

        >>> L = legendre_symbol
        >>> S(45).factors()
        {3: 2, 5: 1}
        >>> jacobi_symbol(7, 45) == L(7, 3)**2 * L(7, 5)**1
        True
    """
    m, n = int_tested(m, n)
    if not n % 2:
        raise ValueError("n should be an odd integer")
    if m < 0 or m > n:
        m = m % n
    if not m:
        return int(n == 1)
    if n == 1 or m == 1:
        return 1
    if igcd(m, n) != 1:
        return 0

    j = 1
    s = trailing(m)
    m = m >> s
    if s % 2 and n % 8 in [3, 5]:
        j *= -1

    while m != 1:
        if m % 4 == 3 and n % 4 == 3:
            j *= -1
        m, n = n % m, m
        s = trailing(m)
        m = m >> s
        if s % 2 and n % 8 in [3, 5]:
            j *= -1
    return j

def totient_(n):
    """returns the number of integers less than n
    and relatively prime to n"""
    if n < 1:
        raise ValueError("n must be a positive integer")
    tot = 0
    for x in xrange(1, n):
        if igcd(x, n) == 1:
            tot += 1
    return tot

def n_order(a,n):
    """ returns the order of a modulo n
    Order of a modulo n is the smallest integer
    k such that a^k leaves a remainder of 1 with n.
    """
    assert igcd(a,n)==1
    if a>n : a=a%n
    for x in xrange(1,totient_(n)+1):
        if (a**x)%n==1:
            return x

def is_primitive_root(a,p):
    """
    returns True if a is a primitive root of p
    """
    assert igcd(a,p) == 1,"The two numbers should be relatively prime"
    if a>p:
        a=a%p
    if n_order(a,p)==totient_(p):
        return True
    else:
        return False

    def get_random_primitive_root(self):
        while True:
            val = random.randint(self._prime // (2 * 2), (self._prime - 1) // 2) * 2 - 1
            if not (val % 3 and val % 5):
                continue

            if igcd(val, self._prime) != 1:
                continue

            if is_primitive_root(val, self._prime):
                return val

def zzx_content(f):
    """Returns integer GCD of coefficients. """
    cont = 0

    for coeff in f:
        cont = igcd(cont, coeff)

        if cont == 1:
            break

    return cont

def is_primitive_root(a, p):
    """
    returns True if a is a primitive root of p
    """
    if igcd(a, p) != 1:
        raise ValueError("The two numbers should be relatively prime")
    if a > p:
        a = a % p
    if n_order(a, p) == totient_(p):
        return True
    else:
        return False

def legendre_symbol(a,p):
    """
    return 1 if a is a quadratic residue of p
    else return -1
    p should be an odd prime by definition
    """
    assert isprime(p) and p!=2,"p should be an odd prime"
    assert igcd(a,p)==1,"The two numbers should be relatively prime"
    if a>p:
        a=a%p
    if is_quad_residue(a,p)==True: return 1
    else : return -1

def _is_nthpow_residue_bign(a, n, m):
    """Returns True if ``x**n == a (mod m)`` has solutions for n > 2."""
    # assert n > 2
    # assert a > 0 and m > 0
    if primitive_root(m) is None:
        # assert m >= 8
        for prime, power in factorint(m).items():
            if not _is_nthpow_residue_bign_prime_power(a, n, prime, power):
                return False
        return True
    f = totient(m)
    k = f // igcd(f, n)
    return pow(a, k, m) == 1

def is_quad_residue(a,p):
    """
    returns True if a is a quadratic residue of p
    p should be a prime and a should be relatively
    prime to p
    """
    assert isprime(p) and p!=2,"p should be an odd prime"
    assert igcd(a,p)==1,"The two numbers should be relatively prime"
    if a>p:
        a=a%p
    rem=(a**((p-1)//2))%p    # a^(p-1 / 2) % p
    if rem==1: return True
    else : return False