def migcdex(x):
     # recursive calcuation of gcd and linear combination
     # for a sequence of integers.
     # Given  (x1, x2, x3)
     # Returns (y1, y1, y3, g)
     # such that g is the gcd and x1*y1+x2*y2+x3*y3 - g = 0
     # Note, that this is only one such linear combination.
     if len(x) == 1:
         return (1, x[0])
     if len(x) == 2:
         return igcdex(x[0], x[-1])
     g = migcdex(x[1:])
     u, v, h = igcdex(x[0], g[-1])
     return tuple([u] + [v*i for i in g[0:-1] ] + [h])

def decipher_elgamal(ct, prk):
    r"""
    Decrypt message with private key

    `ct = (c_{1}, c_{2})`

    `prk = (p, r, d)`

    According to extended Eucliden theorem,
    `u c_{1}^{d} + p n = 1`

    `u \equiv 1/{{c_{1}}^d} \pmod p`

    `u c_{2} \equiv \frac{1}{c_{1}^d} c_{2} \equiv \frac{1}{r^{ad}} c_{2} \pmod p`

    `\frac{1}{r^{ad}} m e^a \equiv \frac{1}{r^{ad}} m {r^{d a}} \equiv m \pmod p`

    Examples
    ========

    >>> from sympy.crypto.crypto import decipher_elgamal
    >>> decipher_elgamal((835, 271), (1031, 14, 636))
    100

    """
    u = igcdex(ct[0] ** prk[2], prk[0])[0]
    return u * ct[1] % prk[0]

def crt(m, v, symmetric=False):
    """Chinese Remainder Theorem.

       The integers in m are assumed to be pairwise coprime.  The output
       is then an integer f, such that f = v_i mod m_i for each pair out
       of v and m.

    """
    mm = 1

    for m_i in m:
        mm *= m_i

    result = 0

    for m_i, v_i in zip(m, v):
        e = mm // m_i
        s, t, g = igcdex(e, m_i)
        c = v_i*s % m_i
        result += c*e

    result %= mm

    if symmetric:
        if result <= mm//2:
            return result
        else:
            return result - mm
    else:
        return result

def _nthroot_mod1(s, q, p, all_roots):
    """
    Root of ``x**q = s mod p``, ``p`` prime and ``q`` divides ``p - 1``

    References
    ==========

    .. [1] A. M. Johnston "A Generalized qth Root Algorithm"

    """
    g = primitive_root(p)
    if not isprime(q):
        r = _nthroot_mod2(s, q, p)
    else:
        f = p - 1
        assert (p - 1) % q == 0
        # determine k
        k = 0
        while f % q == 0:
            k += 1
            f = f // q
        # find z, x, r1
        f1 = igcdex(-f, q)[0] % q
        z = f*f1
        x = (1 + z) // q
        w = pow(g, z, p)
        r1 = pow(s, x, p)
        s1 = pow(s, f, p)
        y = pow(g, f, p)
        h = pow(g, f*q, p)
        t = discrete_log(p, s1, h)
        g2 = pow(g, z*t, p)
        g3 = igcdex(g2, p)[0]
        r = r1*g3 % p
        #assert pow(r, q, p) == s
    res = [r]
    h = pow(g, (p - 1) // q, p)
    #assert pow(h, q, p) == 1
    hx = r
    for i in range(q - 1):
        hx = (hx*h) % p
        res.append(hx)
    if all_roots:
        res.sort()
        return res
    return min(res)

def zp_inv(a, p):
    """Compute multiplicative inverse over GF(p). """
    s, t, g = igcdex(a, p)

    if g == 1:
        return s % p
    else:
        raise ZeroDivisionError("modular division")

def zzx_hensel_lift(p, f, f_list, l):
    """Multifactor Hensel lifting.

       Given a prime p, polynomial f over Z[x] such that lc(f) is a
       unit modulo p,  monic pair-wise coprime polynomials f_i over
       Z[x] satisfying:

                    f = lc(f) f_1 ... f_r (mod p)

       and a positive integer l, returns a list of monic polynomials
       F_1, F_2, ..., F_r satisfying:

                    f = lc(f) F_1 ... F_r (mod p**l)

                    F_i = f_i (mod p), i = 1..r

       For more details on the implemented algorithm refer to:

       [1] J. von zur Gathen, J. Gerhard, Modern Computer Algebra,
           First Edition, Cambridge University Press, 1999, pp. 424

    """
    r = len(f_list)
    lc = zzx_LC(f)

    if r == 1:
        F = zzx_mul_term(f, igcdex(lc, p**l)[0], 0)
        return [ zzx_trunc(F, p**l) ]

    m = p
    k = int(r // 2)
    d = int(ceil(log(l, 2)))

    g = gf_from_int_poly([lc], p)

    for f_i in f_list[:k]:
        g = gf_mul(g, gf_from_int_poly(f_i, p), p)

    h = gf_from_int_poly(f_list[k], p)

    for f_i in f_list[k+1:]:
        h = gf_mul(h, gf_from_int_poly(f_i, p), p)

    s, t, _ = gf_gcdex(g, h, p)

    g = gf_to_int_poly(g, p)
    h = gf_to_int_poly(h, p)
    s = gf_to_int_poly(s, p)
    t = gf_to_int_poly(t, p)

    for _ in range(1, d+1):
        (g, h, s, t), m = zzx_hensel_step(m, f, g, h, s, t), m**2

    return zzx_hensel_lift(p, g, f_list[:k], l) \
         + zzx_hensel_lift(p, h, f_list[k:], l)

def crt1(m):
    """First part of Chinese Remainder Theorem, for multiple application. """
    mm, e, s = 1, [], []

    for m_i in m:
        mm *= m_i

    for m_i in m:
        e.append(mm // m_i)
        s.append(igcdex(e[-1], m_i)[0])

    return mm, e, s

def mod_inverse(a, m):
    """
    Return the number c such that, ( a * c ) % m == 1.
    This means that if a is multiplied by c
    then the remainder obtained on division with m is 1.
    The value of m need not be a prime.
    Return the value of c if some c exists.
    Otherwise, raise an exception stating that the
    modular inverse does not exist.

    Examples
    ========

    >>> from sympy.ntheory.modular import mod_inverse
    >>> mod_inverse(3,11)
    4

    Suppose we wish to find multiplicative inverse x of
    3 modulo 11. This is the same as finding x such
    that 3 * x = 1 (mod 11). One value of x that satisfies
    this congruence is 4. Because 3 * 4 = 12 and 12 = 1 mod(11).

    Similarly

    >>> mod_inverse(5,11)
    9

    Exception Cases

    >>> mod_inverse(2,4)
    Traceback (most recent call last):
    ...
    ValueError: modular inverse does not exist

    References
    ==========
    - https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
    - https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
    """
    a, m = as_int(a), as_int(m)
    x, y, g = igcdex(a, m)
    if g != 1:
        raise ValueError('modular inverse does not exist')
    else:
        return x % m

    def combine(c1, c2):
        """Return the tuple (a, m) which satisfies the requirement
        that n = a + i*m satisfy n = a1 + j*m1 and n = a2 = k*m2.

        Reference:
           http://en.wikipedia.org/wiki/Method_of_successive_substitution
        """
        from sympy.core.numbers import igcdex
        a1, m1 = c1
        a2, m2 = c2
        a, b, c = m1, a2 - a1, m2
        g = reduce(igcd, [a, b, c])
        a, b, c = [i//g for i in [a, b, c]]
        if a != 1:
            inv_a, _, g = igcdex(a, c)
            if g != 1:
                return None
            b *= inv_a
        a, m = a1 + m1*b, m1*c
        return a, m

def test_igcdex():
    assert igcdex(2, 3) == (-1, 1, 1)
    assert igcdex(10, 12) == (-1, 1, 2)
    assert igcdex(100, 2004) == (-20, 1, 4)

def nthroot_mod(a, n, p, all_roots=False):
    """
    find the solutions to ``x**n = a mod p``

    Parameters
    ==========

    a : integer
    n : positive integer
    p : positive integer
    all_roots : if False returns the smallest root, else the list of roots

    Examples
    ========

    >>> from sympy.ntheory.residue_ntheory import nthroot_mod
    >>> nthroot_mod(11, 4, 19)
    8
    >>> nthroot_mod(11, 4, 19, True)
    [8, 11]
    >>> nthroot_mod(68, 3, 109)
    23
    """
    from sympy.core.numbers import igcdex
    if n == 2:
        return sqrt_mod(a, p , all_roots)
    f = totient(p)
    # see Hackman "Elementary Number Theory" (2009), page 76
    if pow(a, f // igcd(f, n), p) != 1:
        return None
    if not isprime(p):
        raise NotImplementedError

    if (p - 1) % n == 0:
        return _nthroot_mod1(a, n, p, all_roots)
    # The roots of ``x**n - a = 0 (mod p)`` are roots of
    # ``gcd(x**n - a, x**(p - 1) - 1) = 0 (mod p)``
    pa = n
    pb = p - 1
    b = 1
    if pa < pb:
        a, pa, b, pb = b, pb, a, pa
    while pb:
        # x**pa - a = 0; x**pb - b = 0
        # x**pa - a = x**(q*pb + r) - a = (x**pb)**q * x**r - a =
        #             b**q * x**r - a; x**r - c = 0; c = b**-q * a mod p
        q, r = divmod(pa, pb)
        c = pow(b, q, p)
        c = igcdex(c, p)[0]
        c = (c * a) % p
        pa, pb = pb, r
        a, b = b, c
    if pa == 1:
        if all_roots:
            res = [a]
        else:
            res = a
    elif pa == 2:
        return sqrt_mod(a, p , all_roots)
    else:
        res = _nthroot_mod1(a, pa, p, all_roots)
    return res

def _sqrt_mod_prime_power(a, p, k):
    """
    find the solutions to ``x**2 = a mod p**k`` when ``a % p != 0``

    Parameters
    ==========

    a : integer
    p : prime number
    k : positive integer

    References
    ==========

    [1] P. Hackman "Elementary Number Theory" (2009),  page 160
    [2] http://www.numbertheory.org/php/squareroot.html
    [3] [Gathen99]_

    Examples
    ========

    >>> from sympy.ntheory.residue_ntheory import _sqrt_mod_prime_power
    >>> _sqrt_mod_prime_power(11, 43, 1)
    [21, 22]
    """
    from sympy.core.numbers import igcdex
    from sympy.polys.domains import ZZ

    pk = p**k
    a = a % pk

    if k == 1:
        if p == 2:
            return [ZZ(a)]
        if not is_quad_residue(a, p):
            return None

        if p % 4 == 3:
            res = pow(a, (p + 1) // 4, p)
        elif p % 8 == 5:
            sign = pow(a, (p - 1) // 4, p)
            if sign == 1:
                res = pow(a, (p + 3) // 8, p)
            else:
                b = pow(4*a, (p - 5) // 8, p)
                x =  (2*a*b) % p
                if pow(x, 2, p) == a:
                    res = x
        else:
            res = _sqrt_mod_tonelli_shanks(a, p)

        # ``_sqrt_mod_tonelli_shanks(a, p)`` is not deterministic;
        # sort to get always the same result
        return sorted([ZZ(res), ZZ(p - res)])

    if k > 1:
        f = factorint(a)
        # see Ref.[2]
        if p == 2:
            if a % 8 != 1:
                return None
            if k <= 3:
               s = set()
               for i in xrange(0, pk, 4):
                    s.add(1 + i)
                    s.add(-1 + i)
               return list(s)
            # according to Ref.[2] for k > 2 there are two solutions
            # (mod 2**k-1), that is four solutions (mod 2**k), which can be
            # obtained from the roots of x**2 = 0 (mod 8)
            rv = [ZZ(1), ZZ(3), ZZ(5), ZZ(7)]
            # hensel lift them to solutions of x**2 = 0 (mod 2**k)
            # if r**2 - a = 0 mod 2**nx but not mod 2**(nx+1)
            # then r + 2**(nx - 1) is a root mod 2**(nx+1)
            n = 3
            res = []
            for r in rv:
                nx = n
                while nx < k:
                    r1 = (r**2 - a) >> nx
                    if r1 % 2:
                        r = r + (1 << (nx - 1))
                    #assert (r**2 - a)% (1 << (nx + 1)) == 0
                    nx += 1
                if r not in res:
                    res.append(r)
                x = r + (1 << (k - 1))
                #assert (x**2 - a) % pk == 0
                if x < (1 << nx) and x not in res:
                    if (x**2 - a) % pk == 0:
                        res.append(x)
            return res
        rv = _sqrt_mod_prime_power(a, p, 1)
        if not rv:
            return None
        r = rv[0]
        fr = r**2 - a
        # hensel lifting with Newton iteration, see Ref.[3] chapter 9
        # with f(x) = x**2 - a; one has f'(a) != 0 (mod p) for p != 2
        n = 1
#.........这里部分代码省略.........