test_numbertheory.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483
  1. import operator
  2. from functools import reduce
  3. import sys
  4. try:
  5. import unittest2 as unittest
  6. except ImportError:
  7. import unittest
  8. import hypothesis.strategies as st
  9. import pytest
  10. from hypothesis import given, settings, example
  11. try:
  12. from hypothesis import HealthCheck
  13. HC_PRESENT = True
  14. except ImportError: # pragma: no cover
  15. HC_PRESENT = False
  16. from .numbertheory import (
  17. SquareRootError,
  18. JacobiError,
  19. factorization,
  20. gcd,
  21. lcm,
  22. jacobi,
  23. inverse_mod,
  24. is_prime,
  25. next_prime,
  26. smallprimes,
  27. square_root_mod_prime,
  28. )
  29. try:
  30. from gmpy2 import mpz
  31. except ImportError:
  32. try:
  33. from gmpy import mpz
  34. except ImportError:
  35. def mpz(x):
  36. return x
  37. BIGPRIMES = (
  38. 999671,
  39. 999683,
  40. 999721,
  41. 999727,
  42. 999749,
  43. 999763,
  44. 999769,
  45. 999773,
  46. 999809,
  47. 999853,
  48. 999863,
  49. 999883,
  50. 999907,
  51. 999917,
  52. 999931,
  53. 999953,
  54. 999959,
  55. 999961,
  56. 999979,
  57. 999983,
  58. )
  59. @pytest.mark.parametrize(
  60. "prime, next_p", [(p, q) for p, q in zip(BIGPRIMES[:-1], BIGPRIMES[1:])]
  61. )
  62. def test_next_prime(prime, next_p):
  63. assert next_prime(prime) == next_p
  64. @pytest.mark.parametrize("val", [-1, 0, 1])
  65. def test_next_prime_with_nums_less_2(val):
  66. assert next_prime(val) == 2
  67. @pytest.mark.slow
  68. @pytest.mark.parametrize("prime", smallprimes)
  69. def test_square_root_mod_prime_for_small_primes(prime):
  70. squares = set()
  71. for num in range(0, 1 + prime // 2):
  72. sq = num * num % prime
  73. squares.add(sq)
  74. root = square_root_mod_prime(sq, prime)
  75. # tested for real with TestNumbertheory.test_square_root_mod_prime
  76. assert root * root % prime == sq
  77. for nonsquare in range(0, prime):
  78. if nonsquare in squares:
  79. continue
  80. with pytest.raises(SquareRootError):
  81. square_root_mod_prime(nonsquare, prime)
  82. def test_square_root_mod_prime_for_2():
  83. a = square_root_mod_prime(1, 2)
  84. assert a == 1
  85. def test_square_root_mod_prime_for_small_prime():
  86. root = square_root_mod_prime(98**2 % 101, 101)
  87. assert root * root % 101 == 9
  88. def test_square_root_mod_prime_for_p_congruent_5():
  89. p = 13
  90. assert p % 8 == 5
  91. root = square_root_mod_prime(3, p)
  92. assert root * root % p == 3
  93. def test_square_root_mod_prime_for_p_congruent_5_large_d():
  94. p = 29
  95. assert p % 8 == 5
  96. root = square_root_mod_prime(4, p)
  97. assert root * root % p == 4
  98. class TestSquareRootModPrime(unittest.TestCase):
  99. def test_power_of_2_p(self):
  100. with self.assertRaises(JacobiError):
  101. square_root_mod_prime(12, 32)
  102. def test_no_square(self):
  103. with self.assertRaises(SquareRootError) as e:
  104. square_root_mod_prime(12, 31)
  105. self.assertIn("no square root", str(e.exception))
  106. def test_non_prime(self):
  107. with self.assertRaises(SquareRootError) as e:
  108. square_root_mod_prime(12, 33)
  109. self.assertIn("p is not prime", str(e.exception))
  110. def test_non_prime_with_negative(self):
  111. with self.assertRaises(SquareRootError) as e:
  112. square_root_mod_prime(697 - 1, 697)
  113. self.assertIn("p is not prime", str(e.exception))
  114. @st.composite
  115. def st_two_nums_rel_prime(draw):
  116. # 521-bit is the biggest curve we operate on, use 1024 for a bit
  117. # of breathing space
  118. mod = draw(st.integers(min_value=2, max_value=2**1024))
  119. num = draw(
  120. st.integers(min_value=1, max_value=mod - 1).filter(
  121. lambda x: gcd(x, mod) == 1
  122. )
  123. )
  124. return num, mod
  125. @st.composite
  126. def st_primes(draw, *args, **kwargs):
  127. if "min_value" not in kwargs: # pragma: no branch
  128. kwargs["min_value"] = 1
  129. prime = draw(
  130. st.sampled_from(smallprimes)
  131. | st.integers(*args, **kwargs).filter(is_prime)
  132. )
  133. return prime
  134. @st.composite
  135. def st_num_square_prime(draw):
  136. prime = draw(st_primes(max_value=2**1024))
  137. num = draw(st.integers(min_value=0, max_value=1 + prime // 2))
  138. sq = num * num % prime
  139. return sq, prime
  140. @st.composite
  141. def st_comp_with_com_fac(draw):
  142. """
  143. Strategy that returns lists of numbers, all having a common factor.
  144. """
  145. primes = draw(
  146. st.lists(st_primes(max_value=2**512), min_size=1, max_size=10)
  147. )
  148. # select random prime(s) that will make the common factor of composites
  149. com_fac_primes = draw(
  150. st.lists(st.sampled_from(primes), min_size=1, max_size=20)
  151. )
  152. com_fac = reduce(operator.mul, com_fac_primes, 1)
  153. # select at most 20 lists (returned numbers),
  154. # each having at most 30 primes (factors) including none (then the number
  155. # will be 1)
  156. comp_primes = draw( # pragma: no branch
  157. st.integers(min_value=1, max_value=20).flatmap(
  158. lambda n: st.lists(
  159. st.lists(st.sampled_from(primes), max_size=30),
  160. min_size=1,
  161. max_size=n,
  162. )
  163. )
  164. )
  165. return [reduce(operator.mul, nums, 1) * com_fac for nums in comp_primes]
  166. @st.composite
  167. def st_comp_no_com_fac(draw):
  168. """
  169. Strategy that returns lists of numbers that don't have a common factor.
  170. """
  171. primes = draw(
  172. st.lists(
  173. st_primes(max_value=2**512), min_size=2, max_size=10, unique=True
  174. )
  175. )
  176. # first select the primes that will create the uncommon factor
  177. # between returned numbers
  178. uncom_fac_primes = draw(
  179. st.lists(
  180. st.sampled_from(primes),
  181. min_size=1,
  182. max_size=len(primes) - 1,
  183. unique=True,
  184. )
  185. )
  186. uncom_fac = reduce(operator.mul, uncom_fac_primes, 1)
  187. # then build composites from leftover primes
  188. leftover_primes = [i for i in primes if i not in uncom_fac_primes]
  189. assert leftover_primes
  190. assert uncom_fac_primes
  191. # select at most 20 lists, each having at most 30 primes
  192. # selected from the leftover_primes list
  193. number_primes = draw( # pragma: no branch
  194. st.integers(min_value=1, max_value=20).flatmap(
  195. lambda n: st.lists(
  196. st.lists(st.sampled_from(leftover_primes), max_size=30),
  197. min_size=1,
  198. max_size=n,
  199. )
  200. )
  201. )
  202. numbers = [reduce(operator.mul, nums, 1) for nums in number_primes]
  203. insert_at = draw(st.integers(min_value=0, max_value=len(numbers)))
  204. numbers.insert(insert_at, uncom_fac)
  205. return numbers
  206. HYP_SETTINGS = {}
  207. if HC_PRESENT: # pragma: no branch
  208. HYP_SETTINGS["suppress_health_check"] = [
  209. HealthCheck.filter_too_much,
  210. HealthCheck.too_slow,
  211. ]
  212. # the factorization() sometimes takes a long time to finish
  213. HYP_SETTINGS["deadline"] = 5000
  214. if "--fast" in sys.argv: # pragma: no cover
  215. HYP_SETTINGS["max_examples"] = 20
  216. HYP_SLOW_SETTINGS = dict(HYP_SETTINGS)
  217. if "--fast" in sys.argv: # pragma: no cover
  218. HYP_SLOW_SETTINGS["max_examples"] = 1
  219. else:
  220. HYP_SLOW_SETTINGS["max_examples"] = 20
  221. class TestIsPrime(unittest.TestCase):
  222. def test_very_small_prime(self):
  223. assert is_prime(23)
  224. def test_very_small_composite(self):
  225. assert not is_prime(22)
  226. def test_small_prime(self):
  227. assert is_prime(123456791)
  228. def test_special_composite(self):
  229. assert not is_prime(10261)
  230. def test_medium_prime_1(self):
  231. # nextPrime[2^256]
  232. assert is_prime(2**256 + 0x129)
  233. def test_medium_prime_2(self):
  234. # nextPrime(2^256+0x129)
  235. assert is_prime(2**256 + 0x12D)
  236. def test_medium_trivial_composite(self):
  237. assert not is_prime(2**256 + 0x130)
  238. def test_medium_non_trivial_composite(self):
  239. assert not is_prime(2**256 + 0x12F)
  240. def test_large_prime(self):
  241. # nextPrime[2^2048]
  242. assert is_prime(mpz(2) ** 2048 + 0x3D5)
  243. def test_pseudoprime_base_19(self):
  244. assert not is_prime(1543267864443420616877677640751301)
  245. def test_pseudoprime_base_300(self):
  246. # F. Arnault "Constructing Carmichael Numbers Which Are Strong
  247. # Pseudoprimes to Several Bases". Journal of Symbolic
  248. # Computation. 20 (2): 151-161. doi:10.1006/jsco.1995.1042.
  249. # Section 4.4 Large Example (a pseudoprime to all bases up to
  250. # 300)
  251. p = int(
  252. "29 674 495 668 685 510 550 154 174 642 905 332 730 "
  253. "771 991 799 853 043 350 995 075 531 276 838 753 171 "
  254. "770 199 594 238 596 428 121 188 033 664 754 218 345 "
  255. "562 493 168 782 883".replace(" ", "")
  256. )
  257. assert is_prime(p)
  258. for _ in range(10):
  259. if not is_prime(p * (313 * (p - 1) + 1) * (353 * (p - 1) + 1)):
  260. break
  261. else:
  262. assert False, "composite not detected"
  263. class TestNumbertheory(unittest.TestCase):
  264. def test_gcd(self):
  265. assert gcd(3 * 5 * 7, 3 * 5 * 11, 3 * 5 * 13) == 3 * 5
  266. assert gcd([3 * 5 * 7, 3 * 5 * 11, 3 * 5 * 13]) == 3 * 5
  267. assert gcd(3) == 3
  268. @unittest.skipUnless(
  269. HC_PRESENT,
  270. "Hypothesis 2.0.0 can't be made tolerant of hard to "
  271. "meet requirements (like `is_prime()`), the test "
  272. "case times-out on it",
  273. )
  274. @settings(**HYP_SLOW_SETTINGS)
  275. @example([877 * 1151, 877 * 1009])
  276. @given(st_comp_with_com_fac())
  277. def test_gcd_with_com_factor(self, numbers):
  278. n = gcd(numbers)
  279. assert 1 in numbers or n != 1
  280. for i in numbers:
  281. assert i % n == 0
  282. @unittest.skipUnless(
  283. HC_PRESENT,
  284. "Hypothesis 2.0.0 can't be made tolerant of hard to "
  285. "meet requirements (like `is_prime()`), the test "
  286. "case times-out on it",
  287. )
  288. @settings(**HYP_SLOW_SETTINGS)
  289. @example([1151, 1069, 1009])
  290. @given(st_comp_no_com_fac())
  291. def test_gcd_with_uncom_factor(self, numbers):
  292. n = gcd(numbers)
  293. assert n == 1
  294. @settings(**HYP_SLOW_SETTINGS)
  295. @given(
  296. st.lists(
  297. st.integers(min_value=1, max_value=2**8192),
  298. min_size=1,
  299. max_size=20,
  300. )
  301. )
  302. def test_gcd_with_random_numbers(self, numbers):
  303. n = gcd(numbers)
  304. for i in numbers:
  305. # check that at least it's a divider
  306. assert i % n == 0
  307. def test_lcm(self):
  308. assert lcm(3, 5 * 3, 7 * 3) == 3 * 5 * 7
  309. assert lcm([3, 5 * 3, 7 * 3]) == 3 * 5 * 7
  310. assert lcm(3) == 3
  311. @settings(**HYP_SLOW_SETTINGS)
  312. @given(
  313. st.lists(
  314. st.integers(min_value=1, max_value=2**8192),
  315. min_size=1,
  316. max_size=20,
  317. )
  318. )
  319. def test_lcm_with_random_numbers(self, numbers):
  320. n = lcm(numbers)
  321. for i in numbers:
  322. assert n % i == 0
  323. @unittest.skipUnless(
  324. HC_PRESENT,
  325. "Hypothesis 2.0.0 can't be made tolerant of hard to "
  326. "meet requirements (like `is_prime()`), the test "
  327. "case times-out on it",
  328. )
  329. @settings(**HYP_SLOW_SETTINGS)
  330. @given(st_num_square_prime())
  331. def test_square_root_mod_prime(self, vals):
  332. square, prime = vals
  333. calc = square_root_mod_prime(square, prime)
  334. assert calc * calc % prime == square
  335. @pytest.mark.slow
  336. @settings(**HYP_SLOW_SETTINGS)
  337. @given(st.integers(min_value=1, max_value=10**12))
  338. @example(265399 * 1526929)
  339. @example(373297**2 * 553991)
  340. def test_factorization(self, num):
  341. factors = factorization(num)
  342. mult = 1
  343. for i in factors:
  344. mult *= i[0] ** i[1]
  345. assert mult == num
  346. def test_factorisation_smallprimes(self):
  347. exp = 101 * 103
  348. assert 101 in smallprimes
  349. assert 103 in smallprimes
  350. factors = factorization(exp)
  351. mult = 1
  352. for i in factors:
  353. mult *= i[0] ** i[1]
  354. assert mult == exp
  355. def test_factorisation_not_smallprimes(self):
  356. exp = 1231 * 1237
  357. assert 1231 not in smallprimes
  358. assert 1237 not in smallprimes
  359. factors = factorization(exp)
  360. mult = 1
  361. for i in factors:
  362. mult *= i[0] ** i[1]
  363. assert mult == exp
  364. def test_jacobi_with_zero(self):
  365. assert jacobi(0, 3) == 0
  366. def test_jacobi_with_one(self):
  367. assert jacobi(1, 3) == 1
  368. @settings(**HYP_SLOW_SETTINGS)
  369. @given(st.integers(min_value=3, max_value=1000).filter(lambda x: x % 2))
  370. def test_jacobi(self, mod):
  371. mod = mpz(mod)
  372. if is_prime(mod):
  373. squares = set()
  374. for root in range(1, mod):
  375. root = mpz(root)
  376. assert jacobi(root * root, mod) == 1
  377. squares.add(root * root % mod)
  378. for i in range(1, mod):
  379. if i not in squares:
  380. i = mpz(i)
  381. assert jacobi(i, mod) == -1
  382. else:
  383. factors = factorization(mod)
  384. for a in range(1, mod):
  385. c = 1
  386. for i in factors:
  387. c *= jacobi(a, i[0]) ** i[1]
  388. assert c == jacobi(a, mod)
  389. @settings(**HYP_SLOW_SETTINGS)
  390. @given(st_two_nums_rel_prime())
  391. def test_inverse_mod(self, nums):
  392. num, mod = nums
  393. inv = inverse_mod(num, mod)
  394. assert 0 < inv < mod
  395. assert num * inv % mod == 1
  396. def test_inverse_mod_with_zero(self):
  397. assert 0 == inverse_mod(0, 11)