test_residue.py 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349
  1. from collections import defaultdict
  2. from sympy.core.containers import Tuple
  3. from sympy.core.singleton import S
  4. from sympy.core.symbol import (Dummy, Symbol)
  5. from sympy.functions.combinatorial.numbers import totient
  6. from sympy.ntheory import n_order, is_primitive_root, is_quad_residue, \
  7. legendre_symbol, jacobi_symbol, primerange, sqrt_mod, \
  8. primitive_root, quadratic_residues, is_nthpow_residue, nthroot_mod, \
  9. sqrt_mod_iter, mobius, discrete_log, quadratic_congruence, \
  10. polynomial_congruence, sieve
  11. from sympy.ntheory.residue_ntheory import _primitive_root_prime_iter, \
  12. _primitive_root_prime_power_iter, _primitive_root_prime_power2_iter, \
  13. _nthroot_mod_prime_power, _discrete_log_trial_mul, _discrete_log_shanks_steps, \
  14. _discrete_log_pollard_rho, _discrete_log_index_calculus, _discrete_log_pohlig_hellman, \
  15. _binomial_mod_prime_power, binomial_mod
  16. from sympy.polys.domains import ZZ
  17. from sympy.testing.pytest import raises
  18. from sympy.core.random import randint, choice
  19. def test_residue():
  20. assert n_order(2, 13) == 12
  21. assert [n_order(a, 7) for a in range(1, 7)] == \
  22. [1, 3, 6, 3, 6, 2]
  23. assert n_order(5, 17) == 16
  24. assert n_order(17, 11) == n_order(6, 11)
  25. assert n_order(101, 119) == 6
  26. assert n_order(11, (10**50 + 151)**2) == 10000000000000000000000000000000000000000000000030100000000000000000000000000000000000000000000022650
  27. raises(ValueError, lambda: n_order(6, 9))
  28. assert is_primitive_root(2, 7) is False
  29. assert is_primitive_root(3, 8) is False
  30. assert is_primitive_root(11, 14) is False
  31. assert is_primitive_root(12, 17) == is_primitive_root(29, 17)
  32. raises(ValueError, lambda: is_primitive_root(3, 6))
  33. for p in primerange(3, 100):
  34. li = list(_primitive_root_prime_iter(p))
  35. assert li[0] == min(li)
  36. for g in li:
  37. assert n_order(g, p) == p - 1
  38. assert len(li) == totient(totient(p))
  39. for e in range(1, 4):
  40. li_power = list(_primitive_root_prime_power_iter(p, e))
  41. li_power2 = list(_primitive_root_prime_power2_iter(p, e))
  42. assert len(li_power) == len(li_power2) == totient(totient(p**e))
  43. assert primitive_root(97) == 5
  44. assert n_order(primitive_root(97, False), 97) == totient(97)
  45. assert primitive_root(97**2) == 5
  46. assert n_order(primitive_root(97**2, False), 97**2) == totient(97**2)
  47. assert primitive_root(40487) == 5
  48. assert n_order(primitive_root(40487, False), 40487) == totient(40487)
  49. # note that primitive_root(40487) + 40487 = 40492 is a primitive root
  50. # of 40487**2, but it is not the smallest
  51. assert primitive_root(40487**2) == 10
  52. assert n_order(primitive_root(40487**2, False), 40487**2) == totient(40487**2)
  53. assert primitive_root(82) == 7
  54. assert n_order(primitive_root(82, False), 82) == totient(82)
  55. p = 10**50 + 151
  56. assert primitive_root(p) == 11
  57. assert n_order(primitive_root(p, False), p) == totient(p)
  58. assert primitive_root(2*p) == 11
  59. assert n_order(primitive_root(2*p, False), 2*p) == totient(2*p)
  60. assert primitive_root(p**2) == 11
  61. assert n_order(primitive_root(p**2, False), p**2) == totient(p**2)
  62. assert primitive_root(4 * 11) is None and primitive_root(4 * 11, False) is None
  63. assert primitive_root(15) is None and primitive_root(15, False) is None
  64. raises(ValueError, lambda: primitive_root(-3))
  65. assert is_quad_residue(3, 7) is False
  66. assert is_quad_residue(10, 13) is True
  67. assert is_quad_residue(12364, 139) == is_quad_residue(12364 % 139, 139)
  68. assert is_quad_residue(207, 251) is True
  69. assert is_quad_residue(0, 1) is True
  70. assert is_quad_residue(1, 1) is True
  71. assert is_quad_residue(0, 2) == is_quad_residue(1, 2) is True
  72. assert is_quad_residue(1, 4) is True
  73. assert is_quad_residue(2, 27) is False
  74. assert is_quad_residue(13122380800, 13604889600) is True
  75. assert [j for j in range(14) if is_quad_residue(j, 14)] == \
  76. [0, 1, 2, 4, 7, 8, 9, 11]
  77. raises(ValueError, lambda: is_quad_residue(1.1, 2))
  78. raises(ValueError, lambda: is_quad_residue(2, 0))
  79. assert quadratic_residues(S.One) == [0]
  80. assert quadratic_residues(1) == [0]
  81. assert quadratic_residues(12) == [0, 1, 4, 9]
  82. assert quadratic_residues(13) == [0, 1, 3, 4, 9, 10, 12]
  83. assert [len(quadratic_residues(i)) for i in range(1, 20)] == \
  84. [1, 2, 2, 2, 3, 4, 4, 3, 4, 6, 6, 4, 7, 8, 6, 4, 9, 8, 10]
  85. assert list(sqrt_mod_iter(6, 2)) == [0]
  86. assert sqrt_mod(3, 13) == 4
  87. assert sqrt_mod(3, -13) == 4
  88. assert sqrt_mod(6, 23) == 11
  89. assert sqrt_mod(345, 690) == 345
  90. assert sqrt_mod(67, 101) == None
  91. assert sqrt_mod(1020, 104729) == None
  92. for p in range(3, 100):
  93. d = defaultdict(list)
  94. for i in range(p):
  95. d[pow(i, 2, p)].append(i)
  96. for i in range(1, p):
  97. it = sqrt_mod_iter(i, p)
  98. v = sqrt_mod(i, p, True)
  99. if v:
  100. v = sorted(v)
  101. assert d[i] == v
  102. else:
  103. assert not d[i]
  104. assert sqrt_mod(9, 27, True) == [3, 6, 12, 15, 21, 24]
  105. assert sqrt_mod(9, 81, True) == [3, 24, 30, 51, 57, 78]
  106. assert sqrt_mod(9, 3**5, True) == [3, 78, 84, 159, 165, 240]
  107. assert sqrt_mod(81, 3**4, True) == [0, 9, 18, 27, 36, 45, 54, 63, 72]
  108. assert sqrt_mod(81, 3**5, True) == [9, 18, 36, 45, 63, 72, 90, 99, 117,\
  109. 126, 144, 153, 171, 180, 198, 207, 225, 234]
  110. assert sqrt_mod(81, 3**6, True) == [9, 72, 90, 153, 171, 234, 252, 315,\
  111. 333, 396, 414, 477, 495, 558, 576, 639, 657, 720]
  112. assert sqrt_mod(81, 3**7, True) == [9, 234, 252, 477, 495, 720, 738, 963,\
  113. 981, 1206, 1224, 1449, 1467, 1692, 1710, 1935, 1953, 2178]
  114. for a, p in [(26214400, 32768000000), (26214400, 16384000000),
  115. (262144, 1048576), (87169610025, 163443018796875),
  116. (22315420166400, 167365651248000000)]:
  117. assert pow(sqrt_mod(a, p), 2, p) == a
  118. n = 70
  119. a, p = 5**2*3**n*2**n, 5**6*3**(n+1)*2**(n+2)
  120. it = sqrt_mod_iter(a, p)
  121. for i in range(10):
  122. assert pow(next(it), 2, p) == a
  123. a, p = 5**2*3**n*2**n, 5**6*3**(n+1)*2**(n+3)
  124. it = sqrt_mod_iter(a, p)
  125. for i in range(2):
  126. assert pow(next(it), 2, p) == a
  127. n = 100
  128. a, p = 5**2*3**n*2**n, 5**6*3**(n+1)*2**(n+1)
  129. it = sqrt_mod_iter(a, p)
  130. for i in range(2):
  131. assert pow(next(it), 2, p) == a
  132. assert type(next(sqrt_mod_iter(9, 27))) is int
  133. assert type(next(sqrt_mod_iter(9, 27, ZZ))) is type(ZZ(1))
  134. assert type(next(sqrt_mod_iter(1, 7, ZZ))) is type(ZZ(1))
  135. assert is_nthpow_residue(2, 1, 5)
  136. #issue 10816
  137. assert is_nthpow_residue(1, 0, 1) is False
  138. assert is_nthpow_residue(1, 0, 2) is True
  139. assert is_nthpow_residue(3, 0, 2) is True
  140. assert is_nthpow_residue(0, 1, 8) is True
  141. assert is_nthpow_residue(2, 3, 2) is True
  142. assert is_nthpow_residue(2, 3, 9) is False
  143. assert is_nthpow_residue(3, 5, 30) is True
  144. assert is_nthpow_residue(21, 11, 20) is True
  145. assert is_nthpow_residue(7, 10, 20) is False
  146. assert is_nthpow_residue(5, 10, 20) is True
  147. assert is_nthpow_residue(3, 10, 48) is False
  148. assert is_nthpow_residue(1, 10, 40) is True
  149. assert is_nthpow_residue(3, 10, 24) is False
  150. assert is_nthpow_residue(1, 10, 24) is True
  151. assert is_nthpow_residue(3, 10, 24) is False
  152. assert is_nthpow_residue(2, 10, 48) is False
  153. assert is_nthpow_residue(81, 3, 972) is False
  154. assert is_nthpow_residue(243, 5, 5103) is True
  155. assert is_nthpow_residue(243, 3, 1240029) is False
  156. assert is_nthpow_residue(36010, 8, 87382) is True
  157. assert is_nthpow_residue(28552, 6, 2218) is True
  158. assert is_nthpow_residue(92712, 9, 50026) is True
  159. x = {pow(i, 56, 1024) for i in range(1024)}
  160. assert {a for a in range(1024) if is_nthpow_residue(a, 56, 1024)} == x
  161. x = { pow(i, 256, 2048) for i in range(2048)}
  162. assert {a for a in range(2048) if is_nthpow_residue(a, 256, 2048)} == x
  163. x = { pow(i, 11, 324000) for i in range(1000)}
  164. assert [ is_nthpow_residue(a, 11, 324000) for a in x]
  165. x = { pow(i, 17, 22217575536) for i in range(1000)}
  166. assert [ is_nthpow_residue(a, 17, 22217575536) for a in x]
  167. assert is_nthpow_residue(676, 3, 5364)
  168. assert is_nthpow_residue(9, 12, 36)
  169. assert is_nthpow_residue(32, 10, 41)
  170. assert is_nthpow_residue(4, 2, 64)
  171. assert is_nthpow_residue(31, 4, 41)
  172. assert not is_nthpow_residue(2, 2, 5)
  173. assert is_nthpow_residue(8547, 12, 10007)
  174. assert is_nthpow_residue(Dummy(even=True) + 3, 3, 2) == True
  175. # _nthroot_mod_prime_power
  176. for p in primerange(2, 10):
  177. for a in range(3):
  178. for n in range(3, 5):
  179. ans = _nthroot_mod_prime_power(a, n, p, 1)
  180. assert isinstance(ans, list)
  181. if len(ans) == 0:
  182. for b in range(p):
  183. assert pow(b, n, p) != a % p
  184. for k in range(2, 10):
  185. assert _nthroot_mod_prime_power(a, n, p, k) == []
  186. else:
  187. for b in range(p):
  188. pred = pow(b, n, p) == a % p
  189. assert not(pred ^ (b in ans))
  190. for k in range(2, 10):
  191. ans = _nthroot_mod_prime_power(a, n, p, k)
  192. if not ans:
  193. break
  194. for b in ans:
  195. assert pow(b, n , p**k) == a
  196. assert nthroot_mod(Dummy(odd=True), 3, 2) == 1
  197. assert nthroot_mod(29, 31, 74) == 45
  198. assert nthroot_mod(1801, 11, 2663) == 44
  199. for a, q, p in [(51922, 2, 203017), (43, 3, 109), (1801, 11, 2663),
  200. (26118163, 1303, 33333347), (1499, 7, 2663), (595, 6, 2663),
  201. (1714, 12, 2663), (28477, 9, 33343)]:
  202. r = nthroot_mod(a, q, p)
  203. assert pow(r, q, p) == a
  204. assert nthroot_mod(11, 3, 109) is None
  205. assert nthroot_mod(16, 5, 36, True) == [4, 22]
  206. assert nthroot_mod(9, 16, 36, True) == [3, 9, 15, 21, 27, 33]
  207. assert nthroot_mod(4, 3, 3249000) is None
  208. assert nthroot_mod(36010, 8, 87382, True) == [40208, 47174]
  209. assert nthroot_mod(0, 12, 37, True) == [0]
  210. assert nthroot_mod(0, 7, 100, True) == [0, 10, 20, 30, 40, 50, 60, 70, 80, 90]
  211. assert nthroot_mod(4, 4, 27, True) == [5, 22]
  212. assert nthroot_mod(4, 4, 121, True) == [19, 102]
  213. assert nthroot_mod(2, 3, 7, True) == []
  214. for p in range(1, 20):
  215. for a in range(p):
  216. for n in range(1, p):
  217. ans = nthroot_mod(a, n, p, True)
  218. assert isinstance(ans, list)
  219. for b in range(p):
  220. pred = pow(b, n, p) == a
  221. assert not(pred ^ (b in ans))
  222. ans2 = nthroot_mod(a, n, p, False)
  223. if ans2 is None:
  224. assert ans == []
  225. else:
  226. assert ans2 in ans
  227. x = Symbol('x', positive=True)
  228. i = Symbol('i', integer=True)
  229. assert _discrete_log_trial_mul(587, 2**7, 2) == 7
  230. assert _discrete_log_trial_mul(941, 7**18, 7) == 18
  231. assert _discrete_log_trial_mul(389, 3**81, 3) == 81
  232. assert _discrete_log_trial_mul(191, 19**123, 19) == 123
  233. assert _discrete_log_shanks_steps(442879, 7**2, 7) == 2
  234. assert _discrete_log_shanks_steps(874323, 5**19, 5) == 19
  235. assert _discrete_log_shanks_steps(6876342, 7**71, 7) == 71
  236. assert _discrete_log_shanks_steps(2456747, 3**321, 3) == 321
  237. assert _discrete_log_pollard_rho(6013199, 2**6, 2, rseed=0) == 6
  238. assert _discrete_log_pollard_rho(6138719, 2**19, 2, rseed=0) == 19
  239. assert _discrete_log_pollard_rho(36721943, 2**40, 2, rseed=0) == 40
  240. assert _discrete_log_pollard_rho(24567899, 3**333, 3, rseed=0) == 333
  241. raises(ValueError, lambda: _discrete_log_pollard_rho(11, 7, 31, rseed=0))
  242. raises(ValueError, lambda: _discrete_log_pollard_rho(227, 3**7, 5, rseed=0))
  243. assert _discrete_log_index_calculus(983, 948, 2, 491) == 183
  244. assert _discrete_log_index_calculus(633383, 21794, 2, 316691) == 68048
  245. assert _discrete_log_index_calculus(941762639, 68822582, 2, 470881319) == 338029275
  246. assert _discrete_log_index_calculus(999231337607, 888188918786, 2, 499615668803) == 142811376514
  247. assert _discrete_log_index_calculus(47747730623, 19410045286, 43425105668, 645239603) == 590504662
  248. assert _discrete_log_pohlig_hellman(98376431, 11**9, 11) == 9
  249. assert _discrete_log_pohlig_hellman(78723213, 11**31, 11) == 31
  250. assert _discrete_log_pohlig_hellman(32942478, 11**98, 11) == 98
  251. assert _discrete_log_pohlig_hellman(14789363, 11**444, 11) == 444
  252. assert discrete_log(1, 0, 2) == 0
  253. raises(ValueError, lambda: discrete_log(-4, 1, 3))
  254. raises(ValueError, lambda: discrete_log(10, 3, 2))
  255. assert discrete_log(587, 2**9, 2) == 9
  256. assert discrete_log(2456747, 3**51, 3) == 51
  257. assert discrete_log(32942478, 11**127, 11) == 127
  258. assert discrete_log(432751500361, 7**324, 7) == 324
  259. assert discrete_log(265390227570863,184500076053622, 2) == 17835221372061
  260. assert discrete_log(22708823198678103974314518195029102158525052496759285596453269189798311427475159776411276642277139650833937,
  261. 17463946429475485293747680247507700244427944625055089103624311227422110546803452417458985046168310373075327,
  262. 123456) == 2068031853682195777930683306640554533145512201725884603914601918777510185469769997054750835368413389728895
  263. args = 5779, 3528, 6215
  264. assert discrete_log(*args) == 687
  265. assert discrete_log(*Tuple(*args)) == 687
  266. assert quadratic_congruence(400, 85, 125, 1600) == [295, 615, 935, 1255, 1575]
  267. assert quadratic_congruence(3, 6, 5, 25) == [3, 20]
  268. assert quadratic_congruence(120, 80, 175, 500) == []
  269. assert quadratic_congruence(15, 14, 7, 2) == [1]
  270. assert quadratic_congruence(8, 15, 7, 29) == [10, 28]
  271. assert quadratic_congruence(160, 200, 300, 461) == [144, 431]
  272. assert quadratic_congruence(100000, 123456, 7415263, 48112959837082048697) == [30417843635344493501, 36001135160550533083]
  273. assert quadratic_congruence(65, 121, 72, 277) == [249, 252]
  274. assert quadratic_congruence(5, 10, 14, 2) == [0]
  275. assert quadratic_congruence(10, 17, 19, 2) == [1]
  276. assert quadratic_congruence(10, 14, 20, 2) == [0, 1]
  277. assert quadratic_congruence(2**48-7, 2**48-1, 4, 2**48) == [8249717183797, 31960993774868]
  278. assert polynomial_congruence(6*x**5 + 10*x**4 + 5*x**3 + x**2 + x + 1,
  279. 972000) == [220999, 242999, 463999, 485999, 706999, 728999, 949999, 971999]
  280. assert polynomial_congruence(x**3 - 10*x**2 + 12*x - 82, 33075) == [30287]
  281. assert polynomial_congruence(x**2 + x + 47, 2401) == [785, 1615]
  282. assert polynomial_congruence(10*x**2 + 14*x + 20, 2) == [0, 1]
  283. assert polynomial_congruence(x**3 + 3, 16) == [5]
  284. assert polynomial_congruence(65*x**2 + 121*x + 72, 277) == [249, 252]
  285. assert polynomial_congruence(x**4 - 4, 27) == [5, 22]
  286. assert polynomial_congruence(35*x**3 - 6*x**2 - 567*x + 2308, 148225) == [86957,
  287. 111157, 122531, 146731]
  288. assert polynomial_congruence(x**16 - 9, 36) == [3, 9, 15, 21, 27, 33]
  289. assert polynomial_congruence(x**6 - 2*x**5 - 35, 6125) == [3257]
  290. raises(ValueError, lambda: polynomial_congruence(x**x, 6125))
  291. raises(ValueError, lambda: polynomial_congruence(x**i, 6125))
  292. raises(ValueError, lambda: polynomial_congruence(0.1*x**2 + 6, 100))
  293. assert binomial_mod(-1, 1, 10) == 0
  294. assert binomial_mod(1, -1, 10) == 0
  295. raises(ValueError, lambda: binomial_mod(2, 1, -1))
  296. assert binomial_mod(51, 10, 10) == 0
  297. assert binomial_mod(10**3, 500, 3**6) == 567
  298. assert binomial_mod(10**18 - 1, 123456789, 4) == 0
  299. assert binomial_mod(10**18, 10**12, (10**5 + 3)**2) == 3744312326
  300. def test_binomial_p_pow():
  301. n, binomials, binomial = 1000, [1], 1
  302. for i in range(1, n + 1):
  303. binomial *= n - i + 1
  304. binomial //= i
  305. binomials.append(binomial)
  306. # Test powers of two, which the algorithm treats slightly differently
  307. trials_2 = 100
  308. for _ in range(trials_2):
  309. m, power = randint(0, n), randint(1, 20)
  310. assert _binomial_mod_prime_power(n, m, 2, power) == binomials[m] % 2**power
  311. # Test against other prime powers
  312. primes = list(sieve.primerange(2*n))
  313. trials = 1000
  314. for _ in range(trials):
  315. m, prime, power = randint(0, n), choice(primes), randint(1, 10)
  316. assert _binomial_mod_prime_power(n, m, prime, power) == binomials[m] % prime**power
  317. def test_deprecated_ntheory_symbolic_functions():
  318. from sympy.testing.pytest import warns_deprecated_sympy
  319. with warns_deprecated_sympy():
  320. assert mobius(3) == -1
  321. with warns_deprecated_sympy():
  322. assert legendre_symbol(2, 3) == -1
  323. with warns_deprecated_sympy():
  324. assert jacobi_symbol(2, 3) == -1