def digits(n, b=10):
    """
    Return a list of the digits of n in base b. The first element in the list
    is b (or -b if n is negative).

    Examples
    ========

    >>> from sympy.ntheory.factor_ import digits
    >>> digits(35)
    [10, 3, 5]
    >>> digits(27, 2)
    [2, 1, 1, 0, 1, 1]
    >>> digits(65536, 256)
    [256, 1, 0, 0]
    >>> digits(-3958, 27)
    [-27, 5, 11, 16]
    """

    b = as_int(b)
    n = as_int(n)
    if b <= 1:
        raise ValueError("b must be >= 2")
    else:
        x, y = abs(n), []
        while x >= b:
            x, r = divmod(x, b)
            y.append(r)
        y.append(x)
        y.append(-b if n < 0 else b)
        y.reverse()
        return y

def is_nthpow_residue(a, n, m):
    """
    Returns True if ``x**n == a (mod m)`` has solutions.

    References
    ==========

    .. [1] P. Hackman "Elementary Number Theory" (2009), page 76

    """
    a, n, m = as_int(a), as_int(n), as_int(m)
    if m <= 0:
        raise ValueError('m must be > 0')
    if n < 0:
        raise ValueError('n must be >= 0')
    if a < 0:
        raise ValueError('a must be >= 0')
    if n == 0:
        if m == 1:
            return False
        return a == 1
    if n == 1:
        return True
    if n == 2:
        return is_quad_residue(a, m)
    return _is_nthpow_residue_bign(a, n, m)

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 legendre_symbol(a, p):
    """
    Returns
    =======

    1. 0 if a is multiple of p
    2. 1 if a is a quadratic residue of p
    3. -1 otherwise

    p should be an odd prime by definition

    Examples
    ========

    >>> from sympy.ntheory import legendre_symbol
    >>> [legendre_symbol(i, 7) for i in range(7)]
    [0, 1, 1, -1, 1, -1, -1]
    >>> list(set([i**2 % 7 for i in range(7)]))
    [0, 1, 2, 4]

    See Also
    ========

    is_quad_residue, jacobi_symbol

    """
    a, p = as_int(a), as_int(p)
    if not isprime(p) or p == 2:
        raise ValueError("p should be an odd prime")
    a = a % p
    if not a:
        return 0
    if is_quad_residue(a, p):
        return 1
    return -1

def is_quad_residue(a, p):
    """
    Returns True if ``a`` (mod ``p``) is in the set of squares mod ``p``,
    i.e a % p in set([i**2 % p for i in range(p)]). If ``p`` is an odd
    prime, an iterative method is used to make the determination:

    >>> from sympy.ntheory import is_quad_residue
    >>> list(set([i**2 % 7 for i in range(7)]))
    [0, 1, 2, 4]
    >>> [j for j in range(7) if is_quad_residue(j, 7)]
    [0, 1, 2, 4]

    See Also
    ========

    legendre_symbol, jacobi_symbol
    """
    a, p = as_int(a), as_int(p)
    if p < 1:
        raise ValueError('p must be > 0')
    if a >= p or a < 0:
        a = a % p
    if a < 2 or p < 3:
        return True
    if not isprime(p):
        if p % 2 and jacobi_symbol(a, p) == -1:
            return False
        r = sqrt_mod(a, p)
        if r is None:
            return False
        else:
            return True

    return pow(a, (p - 1) // 2, p) == 1

def multiplicity(p, n):
    """
    Find the greatest integer m such that p**m divides n.

    Examples
    ========

    >>> from sympy.ntheory import multiplicity
    >>> from sympy.core.numbers import Rational as R
    >>> [multiplicity(5, n) for n in [8, 5, 25, 125, 250]]
    [0, 1, 2, 3, 3]
    >>> multiplicity(3, R(1, 9))
    -2

    """
    try:
            p, n = as_int(p), as_int(n)
    except ValueError:
        if all(isinstance(i, (SYMPY_INTS, Rational)) for i in (p, n)):
            try:
                p = Rational(p)
                n = Rational(n)
                like_min = min(
                    multiplicity(p.p, n.p) if p.p != 1 else oo,
                    multiplicity(p.q, n.q) if p.q != 1 else oo)
                cross_min = min(
                    multiplicity(p.p, n.q) if p.p != 1 else oo,
                    multiplicity(p.q, n.p) if p.q != 1 else oo)
                return like_min - cross_min
            except AttributeError:
                pass
        raise ValueError('expecting ints or fractions, got %s and %s' % (p, n))

    if p == 2:
        return trailing(n)
    if p < 2:
        raise ValueError('p must be an integer, 2 or larger, but got %s' % p)
    if p == n:
        return 1

    m = 0
    n, rem = divmod(n, p)
    while not rem:
        m += 1
        if m > 5:
            # The multiplicity could be very large. Better
            # to increment in powers of two
            e = 2
            while 1:
                ppow = p**e
                if ppow < n:
                    nnew, rem = divmod(n, ppow)
                    if not rem:
                        m += e
                        e *= 2
                        n = nnew
                        continue
                return m + multiplicity(p, n)
        n, rem = divmod(n, p)
    return m

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 = as_int(a), as_int(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.items():
        exponent = group_order
        for f in xrange(e + 1):
            if pow(a, exponent, n) != 1:
                order *= p ** (e - f + 1)
                break
            exponent = exponent // p
    return order

    def search(self, n):
        """Return the indices i, j of the primes that bound n.

        If n is prime then i == j.

        Although n can be an expression, if ceiling cannot convert
        it to an integer then an n error will be raised.

        Examples
        ========

        >>> from sympy import sieve
        >>> sieve.search(25)
        (9, 10)
        >>> sieve.search(23)
        (9, 9)
        """
        from sympy.functions.elementary.integers import ceiling

        # wrapping ceiling in as_int will raise an error if there was a problem
        # determining whether the expression was exactly an integer or not
        test = as_int(ceiling(n))
        n = as_int(n)
        if n < 2:
            raise ValueError("n should be >= 2 but got: %s" % n)
        if n > self._list[-1]:
            self.extend(n)
        b = bisect(self._list, n)
        if self._list[b - 1] == test:
            return b, b
        else:
            return b, b + 1

    def primerange(self, a, b):
        """Generate all prime numbers in the range [a, b).

        Examples
        ========

        >>> from sympy import sieve
        >>> print([i for i in sieve.primerange(7, 18)])
        [7, 11, 13, 17]
        """
        from sympy.functions.elementary.integers import ceiling

        # wrapping ceiling in as_int will raise an error if there was a problem
        # determining whether the expression was exactly an integer or not
        a = max(2, as_int(ceiling(a)))
        b = as_int(ceiling(b))
        if a >= b:
            return
        self.extend(b)
        i = self.search(a)[1]
        maxi = len(self._list) + 1
        while i < maxi:
            p = self._list[i - 1]
            if p < b:
                yield p
                i += 1
            else:
                return

def ones(r, c=None):
    """Returns a matrix of ones with ``r`` rows and ``c`` columns;
    if ``c`` is omitted a square matrix will be returned.

    See Also
    ========

    zeros
    eye
    diag
    """
    from dense import Matrix

    if is_sequence(r):
        SymPyDeprecationWarning(
            feature="The syntax ones([%i, %i])" % tuple(r),
            useinstead="ones(%i, %i)." % tuple(r),
            issue=3381, deprecated_since_version="0.7.2",
        ).warn()
        r, c = r
    else:
        c = r if c is None else c
    r = as_int(r)
    c = as_int(c)
    return Matrix(r, c, [S.One]*r*c)

    def __new__(cls, partition, integer=None):
        """
        Generates a new IntegerPartition object from a list or dictionary.

        The partition can be given as a list of positive integers or a
        dictionary of (integer, multiplicity) items. If the partition is
        preceded by an integer an error will be raised if the partition
        does not sum to that given integer.

        Examples
        ========

        >>> from sympy.combinatorics.partitions import IntegerPartition
        >>> a = IntegerPartition([5, 4, 3, 1, 1])
        >>> a
        IntegerPartition(14, (5, 4, 3, 1, 1))
        >>> print(a)
        [5, 4, 3, 1, 1]
        >>> IntegerPartition({1:3, 2:1})
        IntegerPartition(5, (2, 1, 1, 1))

        If the value that the partition should sum to is given first, a check
        will be made to see n error will be raised if there is a discrepancy:

        >>> IntegerPartition(10, [5, 4, 3, 1])
        Traceback (most recent call last):
        ...
        ValueError: The partition is not valid

        """
        if integer is not None:
            integer, partition = partition, integer
        if isinstance(partition, (dict, Dict)):
            _ = []
            for k, v in sorted(list(partition.items()), reverse=True):
                if not v:
                    continue
                k, v = as_int(k), as_int(v)
                _.extend([k]*v)
            partition = tuple(_)
        else:
            partition = tuple(sorted(map(as_int, partition), reverse=True))
        sum_ok = False
        if integer is None:
            integer = sum(partition)
            sum_ok = True
        else:
            integer = as_int(integer)

        if not sum_ok and sum(partition) != integer:
            raise ValueError("Partition did not add to %s" % integer)
        if any(i < 1 for i in partition):
            raise ValueError("The summands must all be positive.")

        obj = Basic.__new__(cls, integer, partition)
        obj.partition = list(partition)
        obj.integer = integer
        return obj

def sqrt_mod_iter(a, p, domain=int):
    """
    iterate over solutions to ``x**2 = a mod p``

    Parameters
    ==========

    a : integer
    p : positive integer
    domain : integer domain, ``int``, ``ZZ`` or ``Integer``

    Examples
    ========

    >>> from sympy.ntheory.residue_ntheory import sqrt_mod_iter
    >>> list(sqrt_mod_iter(11, 43))
    [21, 22]
    """
    from sympy.polys.galoistools import gf_crt1, gf_crt2
    from sympy.polys.domains import ZZ
    a, p = as_int(a), abs(as_int(p))
    if isprime(p):
        a = a % p
        if a == 0:
            res = _sqrt_mod1(a, p, 1)
        else:
            res = _sqrt_mod_prime_power(a, p, 1)
        if res:
            if domain is ZZ:
                for x in res:
                    yield x
            else:
                for x in res:
                    yield domain(x)
    else:
        f = factorint(p)
        v = []
        pv = []
        for px, ex in f.items():
            if a % px == 0:
                rx = _sqrt_mod1(a, px, ex)
                if not rx:
                    raise StopIteration
            else:
                rx = _sqrt_mod_prime_power(a, px, ex)
                if not rx:
                    raise StopIteration
            v.append(rx)
            pv.append(px**ex)
        mm, e, s = gf_crt1(pv, ZZ)
        if domain is ZZ:
            for vx in _product(*v):
                r = gf_crt2(vx, pv, mm, e, s, ZZ)
                yield r
        else:
            for vx in _product(*v):
                r = gf_crt2(vx, pv, mm, e, s, ZZ)
                yield domain(r)

def _number_theoretic_transform(seq, prime, inverse=False):
    """Utility function for the Number Theoretic Transform"""

    if not iterable(seq):
        raise TypeError("Expected a sequence of integer coefficients "
                        "for Number Theoretic Transform")

    p = as_int(prime)
    if isprime(p) == False:
        raise ValueError("Expected prime modulus for "
                        "Number Theoretic Transform")

    a = [as_int(x) % p for x in seq]

    n = len(a)
    if n < 1:
        return a

    b = n.bit_length() - 1
    if n&(n - 1):
        b += 1
        n = 2**b

    if (p - 1) % n:
        raise ValueError("Expected prime modulus of the form (m*2**k + 1)")

    a += [0]*(n - len(a))
    for i in range(1, n):
        j = int(ibin(i, b, str=True)[::-1], 2)
        if i < j:
            a[i], a[j] = a[j], a[i]

    pr = primitive_root(p)

    rt = pow(pr, (p - 1) // n, p)
    if inverse:
        rt = pow(rt, p - 2, p)

    w = [1]*(n // 2)
    for i in range(1, n // 2):
        w[i] = w[i - 1]*rt % p

    h = 2
    while h <= n:
        hf, ut = h // 2, n // h
        for i in range(0, n, h):
            for j in range(hf):
                u, v = a[i + j], a[i + j + hf]*w[ut * j]
                a[i + j], a[i + j + hf] = (u + v) % p, (u - v) % p
        h *= 2

    if inverse:
        rv = pow(n, p - 2, p)
        a = [x*rv % p for x in a]

    return a

 def as_dict(e):
     # creates a dictionary of the expression using either
     # as_coefficients_dict or as_powers_dict, depending on Func
     # SHARES meth
     d = getattr(e, meth, lambda: {a: S.One for a in e.args})()
     for k in list(d.keys()):
         try:
             as_int(d[k])
         except ValueError:
             d[F(k, d.pop(k))] = S.One
     return d

 def give(a, b, seq=seed):
     a, b = as_int(a), as_int(b)
     w = b - a
     if w < 0:
         raise ValueError('_randint got empty range')
     try:
         x = seq.pop()
     except IndexError:
         raise ValueError('_randint sequence was too short')
     if a <= x <= b:
         return x
     else:
         return give(a, b, seq)

    def __new__(cls, *args, **kwargs):
        self = object.__new__(cls)
        if len(args) == 1 and isinstance(args[0], SparseMatrix):
            self.rows = args[0].rows
            self.cols = args[0].cols
            self._smat = dict(args[0]._smat)
            return self

        self._smat = {}

        if len(args) == 3:
            self.rows = as_int(args[0])
            self.cols = as_int(args[1])

            if isinstance(args[2], Callable):
                op = args[2]
                for i in range(self.rows):
                    for j in range(self.cols):
                        value = self._sympify(
                            op(self._sympify(i), self._sympify(j)))
                        if value:
                            self._smat[(i, j)] = value
            elif isinstance(args[2], (dict, Dict)):
                # manual copy, copy.deepcopy() doesn't work
                for key in args[2].keys():
                    v = args[2][key]
                    if v:
                        self._smat[key] = self._sympify(v)
            elif is_sequence(args[2]):
                if len(args[2]) != self.rows*self.cols:
                    raise ValueError(
                        'List length (%s) != rows*columns (%s)' %
                        (len(args[2]), self.rows*self.cols))
                flat_list = args[2]
                for i in range(self.rows):
                    for j in range(self.cols):
                        value = self._sympify(flat_list[i*self.cols + j])
                        if value:
                            self._smat[(i, j)] = value
        else:
            # handle full matrix forms with _handle_creation_inputs
            r, c, _list = Matrix._handle_creation_inputs(*args)
            self.rows = r
            self.cols = c
            for i in range(self.rows):
                for j in range(self.cols):
                    value = _list[self.cols*i + j]
                    if value:
                        self._smat[(i, j)] = value
        return self

 def zeros(cls, r, c=None):
     """Return an r x c matrix of zeros, square if c is omitted."""
     if is_sequence(r):
         SymPyDeprecationWarning(
             feature="The syntax zeros([%i, %i])" % tuple(r),
             useinstead="zeros(%i, %i)." % tuple(r),
             issue=3381, deprecated_since_version="0.7.2",
         ).warn()
         r, c = r
     else:
         c = r if c is None else c
     r = as_int(r)
     c = as_int(c)
     return cls._new(r, c, [S.Zero]*r*c)

    def __new__(cls, *args):
        from sympy.functions.elementary.integers import ceiling
        # expand range
        slc = slice(*args)
        start, stop, step = slc.start or 0, slc.stop, slc.step or 1
        try:
            start, stop, step = [w if w in [S.NegativeInfinity, S.Infinity] else S(as_int(w))
                                 for w in (start, stop, step)]
        except ValueError:
            raise ValueError("Inputs to Range must be Integer Valued\n" +
                    "Use ImageSets of Ranges for other cases")

        if not step.is_finite:
            raise ValueError("Infinite step is not allowed")
        if start == stop:
            return S.EmptySet

        n = ceiling((stop - start)/step)
        if n <= 0:
            return S.EmptySet

        # normalize args: regardless of how they are entered they will show
        # canonically as Range(inf, sup, step) with step > 0
        if n.is_finite:
            start, stop = sorted((start, start + (n - 1)*step))
        else:
            start, stop = sorted((start, stop - step))

        step = abs(step)
        if (start, stop) == (S.NegativeInfinity, S.Infinity):
            raise ValueError("Both the start and end value of "
                             "Range cannot be unbounded")
        else:
            return Basic.__new__(cls, start, stop + step, step)

def _primitive_root_prime_iter(p):
    """
    Generates the primitive roots for a prime ``p``

    References
    ==========

    [1] W. Stein "Elementary Number Theory" (2011), page 44

    Examples
    ========

    >>> from sympy.ntheory.residue_ntheory import _primitive_root_prime_iter
    >>> list(_primitive_root_prime_iter(19))
    [2, 3, 10, 13, 14, 15]
    """
    p = as_int(p)
    v = [(p - 1) // i for i in factorint(p - 1).keys()]
    a = 2
    while a < p:
        for pw in v:
            if pow(a, pw, p) == 1:
                break
        else:
            yield a
        a += 1

def antidivisor_count(n):
    """
    Return the number of antidivisors [1]_ of ``n``.

    References
    ==========

    .. [1] formula from https://oeis.org/A066272

    Examples
    ========

    >>> from sympy.ntheory.factor_ import antidivisor_count
    >>> antidivisor_count(13)
    4
    >>> antidivisor_count(27)
    5

    See Also
    ========

    factorint, divisors, antidivisors, divisor_count, totient
    """

    n = as_int(abs(n))
    if n <= 2:
        return 0
    return divisor_count(2*n-1) + divisor_count(2*n+1) + \
        divisor_count(n) - divisor_count(n, 2) - 5