def roots_rational(f):
    """Returns a list of rational roots of a polynomial."""

    try:
        g = f.as_integer()[1]
    except CoefficientError:
        return []

    LC_divs = divisors(int(g.LC))
    TC_divs = divisors(int(g.TC))

    if not g(S.Zero):
        zeros = [S.Zero]
    else:
        zeros = []

    for p in LC_divs:
        for q in TC_divs:
            zero = Rational(p, q)

            if not g(zero):
                zeros.append(zero)

            if not g(-zero):
                zeros.append(-zero)

    return zeros

def roots_rational(f):
    """Returns a list of rational roots of a polynomial."""
    domain = f.get_domain()

    if domain.is_QQ:
        _, f = f.clear_denoms()
    elif domain.is_ZZ:
        f = f.set_domain('QQ')
    else:
        return []

    LC_divs = divisors(int(f.LC()))
    EC_divs = divisors(int(f.EC()))

    if not f.eval(S.Zero):
        zeros = [S.Zero]
    else:
        zeros = []

    for p in LC_divs:
        for q in EC_divs:
            zero = Rational(p, q)

            if not f.eval(zero):
                zeros.append(zero)

            if not f.eval(-zero):
                zeros.append(-zero)

    return sorted(zeros, key=default_sort_key)

def _integer_basis(poly):
    """Compute coefficient basis for a polynomial over integers. """
    monoms, coeffs = zip(*poly.terms())

    monoms, = zip(*monoms)
    coeffs = map(abs, coeffs)

    if coeffs[0] < coeffs[-1]:
        coeffs = list(reversed(coeffs))
    else:
        return None

    monoms = monoms[:-1]
    coeffs = coeffs[:-1]

    divs = reversed(divisors(gcd_list(coeffs))[1:])

    try:
        div = divs.next()
    except StopIteration:
        return None

    while True:
        for monom, coeff in zip(monoms, coeffs):
            if coeff % div**monom != 0:
                try:
                    div = divs.next()
                except StopIteration:
                    return None
                else:
                    break
        else:
            return div

def calc_recurrence_seq(p):
    decimal.getcontext().prec = p - 1
    r = reciprocal(p)
    for d in ntheory.divisors(p - 1):
        if d == 1 or d == p - 1:
            continue
        if r[:d] == r[d:2 * d]:
            return r[:d]
    return r[:p - 1]

def MultiOrder(n):
    div = divisors(n)
    if list(set(div).intersection((2,5,10))) != []:
        return None
    order=1
    product = 10
    while product != 1:
        order+=1
        product = (product*10)%n
    return order

def _integer_basis(poly):
    """Compute coefficient basis for a polynomial over integers.

    Returns the integer ``div`` such that substituting ``x = div*y``
    ``p(x) = m*q(y)`` where the coefficients of ``q`` are smaller
    than those of ``p``.

    For example ``x**5 + 512*x + 1024 = 0``
    with ``div = 4`` becomes ``y**5 + 2*y + 1 = 0``

    Returns the integer ``div`` or ``None`` if there is no possible scaling.

    Examples
    ========

    >>> from sympy.polys import Poly
    >>> from sympy.abc import x
    >>> from sympy.polys.polyroots import _integer_basis
    >>> p = Poly(x**5 + 512*x + 1024, x, domain='ZZ')
    >>> _integer_basis(p)
    4
    """
    monoms, coeffs = list(zip(*poly.terms()))

    monoms, = list(zip(*monoms))
    coeffs = list(map(abs, coeffs))

    if coeffs[0] < coeffs[-1]:
        coeffs = list(reversed(coeffs))
        n = monoms[0]
        monoms = [n - i for i in reversed(monoms)]
    else:
        return None

    monoms = monoms[:-1]
    coeffs = coeffs[:-1]

    divs = reversed(divisors(gcd_list(coeffs))[1:])

    try:
        div = next(divs)
    except StopIteration:
        return None

    while True:
        for monom, coeff in zip(monoms, coeffs):
            if coeff % div**monom != 0:
                try:
                    div = next(divs)
                except StopIteration:
                    return None
                else:
                    break
        else:
            return div

def divisors_and_divisor_count():
    assert divisors(-1) == [1]
    assert divisors(0) == []
    assert divisors(1) == [1]
    assert divisors(2) == [1, 2]
    assert divisors(3) == [1, 3]
    assert divisors(17) == [1, 17]
    assert divisors(10) == [1, 2, 5, 10]
    assert divisors(100) == [1, 2, 4, 5, 10, 20, 25, 50, 100]
    assert divisors(101) == [1, 101]

    assert divisor_count(0) == 0
    assert divisor_count(-1) == 1
    assert divisor_count(1) == 1
    assert divisor_count(6) == 4
    assert divisor_count(12) == 6

def find_rational_roots(f):
    if not all(x.is_rational for x in f.all_coeffs()):
        return []
    a = f.nth(0)
    b = f.nth(f.degree())
    result = []
    for p in divisors(a):
        for q in divisors(b):
            if igcd(p, q) == 1:
                if f(Rational(p, q)) == 0:
                    result += [Rational(p, q)]
                if f(Rational(-p, q)) == 0:
                    result += [Rational(-p, q)]
    if len(result) > 0:
        add_comment("Use the rational root test")
        add_comment("Test all rational numbers in the form p/q")
        add_comment("where p is a divisor of {} and q is a divisor of {}", str(a), str(b))
        add_comment("We have the following roots")
        for r in result:
            add_eq(f.gen, r)
    return result

def test_divisors_and_divisor_count():
    assert divisors(-1) == [1]
    assert divisors(0) == []
    assert divisors(1) == [1]
    assert divisors(2) == [1, 2]
    assert divisors(3) == [1, 3]
    assert divisors(17) == [1, 17]
    assert divisors(10) == [1, 2, 5, 10]
    assert divisors(100) == [1, 2, 4, 5, 10, 20, 25, 50, 100]
    assert divisors(101) == [1, 101]

    assert divisor_count(0) == 0
    assert divisor_count(-1) == 1
    assert divisor_count(1) == 1
    assert divisor_count(6) == 4
    assert divisor_count(12) == 6

    assert divisor_count(180, 3) == divisor_count(180//3)
    assert divisor_count(2*3*5, 7) == 0

def _inv_totient_estimate(m):
    """
    Find ``(L, U)`` such that ``L <= phi^-1(m) <= U``.

    Examples
    ========

    >>> from sympy.polys.polyroots import _inv_totient_estimate

    >>> _inv_totient_estimate(192)
    (192, 840)
    >>> _inv_totient_estimate(400)
    (400, 1750)

    """
    primes = [ d + 1 for d in divisors(m) if isprime(d + 1) ]

    a, b = 1, 1

    for p in primes:
        a *= p
        b *= p - 1

    L = m
    U = int(math.ceil(m*(float(a)/b)))

    P = p = 2
    primes = []

    while P <= U:
        p = nextprime(p)
        primes.append(p)
        P *= p

    P //= p
    b = 1

    for p in primes[:-1]:
        b *= p - 1

    U = int(math.ceil(m*(float(P)/b)))

    return L, U

def test_factorint():
    assert sorted(factorint(123456).items()) == [(2, 6), (3, 1), (643, 1)]
    assert primefactors(123456) == [2, 3, 643]
    assert factorint(-16) == {-1:1, 2:4}
    assert factorint(2**(2**6) + 1) == {274177:1, 67280421310721:1}
    assert factorint(5951757) == {3:1, 7:1, 29:2, 337:1}
    assert factorint(64015937) == {7993:1, 8009:1}
    assert divisors(1) == [1]
    assert divisors(2) == [1, 2]
    assert divisors(3) == [1, 3]
    assert divisors(10) == [1, 2, 5, 10]
    assert divisors(100) == [1, 2, 4, 5, 10, 20, 25, 50, 100]
    assert divisors(101) == [1, 101]
    assert factorint(0) == {0:1}
    assert factorint(1) == {}
    assert factorint(-1) == {-1:1}
    assert factorint(-2) == {-1:1, 2:1}
    assert factorint(-16) == {-1:1, 2:4}
    assert factorint(2) == {2:1}
    assert factorint(126) == {2:1, 3:2, 7:1}
    assert factorint(123456) == {2:6, 3:1, 643:1}
    assert factorint(5951757) == {3:1, 7:1, 29:2, 337:1}
    assert factorint(64015937) == {7993:1, 8009:1}
    assert factorint(2**(2**6) + 1) == {274177:1, 67280421310721:1}
    assert multiproduct(factorint(fac(200))) == fac(200)
    for b, e in factorint(fac(150)).items():
        assert e == fac_multiplicity(150, b)
    assert factorint(103005006059**7) == {103005006059:7}
    assert factorint(31337**191) == {31337:191}
    assert factorint(2**1000 * 3**500 * 257**127 * 383**60) == \
        {2:1000, 3:500, 257:127, 383:60}
    assert len(factorint(fac(10000))) == 1229
    assert factorint(12932983746293756928584532764589230) == \
        {2: 1, 5: 1, 73: 1, 727719592270351: 1, 63564265087747: 1, 383: 1}
    assert factorint(727719592270351) == {727719592270351:1}
    assert factorint(2**64+1, use_trial=False) == factorint(2**64+1)
    for n in range(60000):
        assert multiproduct(factorint(n)) == n
    assert pollard_rho(2**64+1, seed=1) == 274177
    assert pollard_rho(19, seed=1) is None
    assert factorint(3,2) == {3: 1}
    assert factorint(12345) == {3: 1, 5: 1, 823: 1}
    assert factorint(12345, 3) == {12345: 1} # there are no factors less than 3

def test_factor():
    assert trial(1) == []
    assert trial(2) == [(2,1)]
    assert trial(3) == [(3,1)]
    assert trial(4) == [(2,2)]
    assert trial(5) == [(5,1)]
    assert trial(128) == [(2,7)]
    assert trial(720) == [(2,4), (3,2), (5,1)]
    assert factorint(123456) == [(2, 6), (3, 1), (643, 1)]
    assert primefactors(123456) == [2, 3, 643]
    assert factorint(-16) == [(-1, 1), (2, 4)]
    assert factorint(2**(2**6) + 1) == [(274177, 1), (67280421310721, 1)]
    assert factorint(5951757) == [(3, 1), (7, 1), (29, 2), (337, 1)]
    assert factorint(64015937) == [(7993, 1), (8009, 1)]
    assert divisors(1) == [1]
    assert divisors(2) == [1, 2]
    assert divisors(3) == [1, 3]
    assert divisors(10) == [1, 2, 5, 10]
    assert divisors(100) == [1, 2, 4, 5, 10, 20, 25, 50, 100]
    assert divisors(101) == [1, 101]
    assert pollard_rho(2**64+1, max_iters=1, seed=1) == 274177
    assert pollard_rho(19) is None

def test_divisors():
    assert divisors(28) == [1, 2, 4, 7, 14, 28]
    assert [x for x in divisors(3*5*7, 1)] == [1, 3, 5, 15, 7, 21, 35, 105]
    assert divisors(0) == []

def test_issue_6981():
    S = set(divisors(4)).union(set(divisors(Integer(2))))
    assert S == set([1,2,4])

from sympy.ntheory import divisors
import numpy as np


sum = 0

for i in xrange(1,10000):
    divSum = np.sum(divisors(i)[:-1])
    if i == np.sum(divisors(divSum)[:-1]) and i !=divSum:
        sum += i

print sum

def calc_some(p):
    if p <= 5:
        raise ValueError("p must be larger than 5")
    if not sympy.isprime(p):
        raise ValueError("specify prime. actual p = {}, ntheory.divisors = {}".format(p, ntheory.divisors(p)))
    r = reciprocal(p)
    seq = calc_recurrence_seq(p)
    d = len(seq)
    former = r[:d // 2]
    latter = r[d // 2:d]
    if d % 2 == 0:
        print('length of recurrence is even')
        print('1/{} = '.format(p), r)
        print('former = ', former)
        print('latter = ', latter)
        print('former + latter =\n', Decimal(former) + Decimal(latter))
    else:
        print('length of recurrence is not even. In fact')
        print('1/{}'.format(p), r)
        print('seq = ', seq)
        print('len(seq) = ', len(seq))
        print('But you may find something')
        if latter[0] == 9:
            print("you're lucky")
            print('former = ', former)
            print('latter = ', latter)
            print('former + latter =\n', Decimal(former) + Decimal(latter))