Source code for sympy.ntheory.residue_ntheory
from __future__ import print_function, division
from sympy.core.numbers import igcd, igcdex
from sympy.core.compatibility import as_int, xrange
from .primetest import isprime
from .factor_ import factorint, trailing, totient
from random import randint
[docs]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
"""
from collections import defaultdict
a, n = as_int(a), as_int(n)
if igcd(a, n) != 1:
raise ValueError("The two numbers should be relatively prime")
factors = defaultdict(int)
f = factorint(n)
for px, kx in f.items():
if kx > 1:
factors[px] += kx - 1
fpx = factorint(px - 1)
for py, ky in fpx.items():
factors[py] += ky
group_order = 1
for px, kx in factors.items():
group_order *= px**kx
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 _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
[docs]def primitive_root(p):
"""
Returns the smallest primitive root or None
References
==========
[1] W. Stein "Elementary Number Theory" (2011), page 44
[2] P. Hackman "Elementary Number Theory" (2009), Chapter C
Parameters
==========
p : positive integer
Examples
========
>>> from sympy.ntheory.residue_ntheory import primitive_root
>>> primitive_root(19)
2
"""
p = as_int(p)
if p < 1:
raise ValueError('p is required to be positive')
if p <= 2:
return 1
f = factorint(p)
if len(f) > 2:
return None
if len(f) == 2:
if 2 not in f or f[2] > 1:
return None
# case p = 2*p1**k, p1 prime
for p1, e1 in f.items():
if p1 != 2:
break
i = 1
while i < p:
i += 2
if i % p1 == 0:
continue
if is_primitive_root(i, p):
return i
else:
if 2 in f:
if p == 4:
return 3
return None
p1, n = list(f.items())[0]
if n > 1:
# see Ref [2], page 81
g = primitive_root(p1)
if is_primitive_root(g, p1**2):
return g
else:
for i in xrange(2, g + p1 + 1):
if igcd(i, p) == 1 and is_primitive_root(i, p):
return i
return next(_primitive_root_prime_iter(p))
[docs]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 _sqrt_mod_tonelli_shanks(a, p):
"""
Returns the square root in the case of ``p`` prime with ``p == 1 (mod 8)``
References
==========
R. Crandall and C. Pomerance "Prime Numbers", 2nt Ed., page 101
"""
s = trailing(p - 1)
t = p >> s
# find a non-quadratic residue
while 1:
d = randint(2, p - 1)
r = legendre_symbol(d, p)
if r == -1:
break
#assert legendre_symbol(d, p) == -1
A = pow(a, t, p)
D = pow(d, t, p)
m = 0
for i in xrange(s):
adm = A*pow(D, m, p) % p
adm = pow(adm, 2**(s - 1 - i), p)
if adm % p == p - 1:
m += 2**i
#assert A*pow(D, m, p) % p == 1
x = pow(a, (t + 1)//2, p)*pow(D, m//2, p) % p
return x
[docs]def sqrt_mod(a, p, all_roots=False):
"""
find a root of ``x**2 = a mod p``
Parameters
==========
a : integer
p : positive integer
all_roots : if True the list of roots is returned or None
Notes
=====
If there is no root it is returned None; else the returned root
is less or equal to ``p // 2``; in general is not the smallest one.
It is returned ``p // 2`` only if it is the only root.
Use ``all_roots`` only when it is expected that all the roots fit
in memory; otherwise use ``sqrt_mod_iter``.
Examples
========
>>> from sympy.ntheory import sqrt_mod
>>> sqrt_mod(11, 43)
21
>>> sqrt_mod(17, 32, True)
[7, 9, 23, 25]
"""
if all_roots:
return sorted(list(sqrt_mod_iter(a, p)))
try:
p = abs(as_int(p))
it = sqrt_mod_iter(a, p)
r = next(it)
if r > p // 2:
return p - r
elif r < p // 2:
return r
else:
try:
r = next(it)
if r > p // 2:
return p - r
except StopIteration:
pass
return r
except StopIteration:
return None
def _product(*iters):
"""
cartesian product generator
Notes
=====
Unlike itertools.product, it works also with iterables which do not fit
in memory. See http://bugs.python.org/issue10109
Author: Fernando Sumudu
with small changes
"""
import itertools
inf_iters = tuple(itertools.cycle(enumerate(it)) for it in iters)
num_iters = len(inf_iters)
cur_val = [None]*num_iters
first_v = True
while True:
i, p = 0, num_iters
while p and not i:
p -= 1
i, cur_val[p] = next(inf_iters[p])
if not p and not i:
if first_v:
first_v = False
else:
break
yield cur_val
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_crt, 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 _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
px = p
while 1:
n1 = n
n1 *= 2
if n1 > k:
break
n = n1
px = px**2
frinv = igcdex(2*r, px)[0]
r = (r - fr*frinv) % px
fr = r**2 - a
if n < k:
px = p**k
frinv = igcdex(2*r, px)[0]
r = (r - fr*frinv) % px
return [r, px - r]
def _sqrt_mod1(a, p, n):
"""
find solution to ``x**2 == a mod p**n`` when ``a % p == 0``
see http://www.numbertheory.org/php/squareroot.html
"""
pn = p**n
a = a % pn
if a == 0:
# case gcd(a, p**k) = p**n
m = n // 2
if n % 2 == 1:
pm1 = p**(m + 1)
def _iter0a():
i = 0
while i < pn:
yield i
i += pm1
return _iter0a()
else:
pm = p**m
def _iter0b():
i = 0
while i < pn:
yield i
i += pm
return _iter0b()
# case gcd(a, p**k) = p**r, r < n
f = factorint(a)
r = f[p]
if r % 2 == 1:
return None
m = r // 2
a1 = a >> r
if p == 2:
if n - r == 1:
pnm1 = 1 << (n - m + 1)
pm1 = 1 << (m + 1)
def _iter1():
k = 1 << (m + 2)
i = 1 << m
while i < pnm1:
j = i
while j < pn:
yield j
j += k
i += pm1
return _iter1()
if n - r == 2:
res = _sqrt_mod_prime_power(a1, p, n - r)
if res is None:
return None
pnm = 1 << (n - m)
def _iter2():
s = set()
for r in res:
i = 0
while i < pn:
x = (r << m) + i
if x not in s:
s.add(x)
yield x
i += pnm
return _iter2()
if n - r > 2:
res = _sqrt_mod_prime_power(a1, p, n - r)
if res is None:
return None
pnm1 = 1 << (n - m - 1)
def _iter3():
s = set()
for r in res:
i = 0
while i < pn:
x = ((r << m) + i) % pn
if x not in s:
s.add(x)
yield x
i += pnm1
return _iter3()
else:
m = r // 2
a1 = a // p**r
res1 = _sqrt_mod_prime_power(a1, p, n - r)
if res1 is None:
return None
pm = p**m
pnr = p**(n-r)
pnm = p**(n-m)
def _iter4():
s = set()
pm = p**m
for rx in res1:
i = 0
while i < pnm:
x = ((rx + i) % pn)
if x not in s:
s.add(x)
yield x*pm
i += pnr
return _iter4()
[docs]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
[docs]def is_nthpow_residue(a, n, m):
"""
Returns True if ``x**n == a (mod m)`` has solutions.
References
==========
P. Hackman "Elementary Number Theory" (2009), page 76
"""
if n == 1:
return True
if n == 2:
return is_quad_residue(a, m)
f = totient(m)
k = f // igcd(f, n)
return pow(a, k, m) == 1
def _nthroot_mod2(s, q, p):
f = factorint(q)
v = []
for b, e in f.items():
v.extend([b]*e)
for qx in v:
s = _nthroot_mod1(s, qx, p, False)
return s
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)
# find t discrete log of s1 base h, h**x = s1 mod p
# used a naive implementation
# TODO implement using Ref [1]
pr = 1
for t in xrange(p):
if pr == s1:
break
pr = pr*h % p
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)
[docs]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
[docs]def quadratic_residues(p):
"""
Returns the list of quadratic residues.
Examples
========
>>> from sympy.ntheory.residue_ntheory import quadratic_residues
>>> quadratic_residues(7)
[0, 1, 2, 4]
"""
r = set()
for i in xrange(p // 2 + 1):
r.add(pow(i, 2, p))
return sorted(list(r))
[docs]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
[docs]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
See Also
========
is_quad_residue, legendre_symbol
"""
m, n = as_int(m), as_int(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
