test_primetest.py 9.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235
  1. from math import gcd
  2. from sympy.ntheory.generate import Sieve, sieve
  3. from sympy.ntheory.primetest import (mr, _lucas_extrastrong_params, is_lucas_prp, is_square,
  4. is_strong_lucas_prp, is_extra_strong_lucas_prp,
  5. proth_test, isprime, is_euler_pseudoprime,
  6. is_gaussian_prime, is_fermat_pseudoprime, is_euler_jacobi_pseudoprime,
  7. MERSENNE_PRIME_EXPONENTS, _lucas_lehmer_primality_test,
  8. is_mersenne_prime)
  9. from sympy.testing.pytest import slow, raises
  10. from sympy.core.numbers import I, Float
  11. def test_is_fermat_pseudoprime():
  12. assert is_fermat_pseudoprime(5, 1)
  13. assert is_fermat_pseudoprime(9, 1)
  14. def test_euler_pseudoprimes():
  15. assert is_euler_pseudoprime(13, 1)
  16. assert is_euler_pseudoprime(15, 1)
  17. assert is_euler_pseudoprime(17, 6)
  18. assert is_euler_pseudoprime(101, 7)
  19. assert is_euler_pseudoprime(1009, 10)
  20. assert is_euler_pseudoprime(11287, 41)
  21. raises(ValueError, lambda: is_euler_pseudoprime(0, 4))
  22. raises(ValueError, lambda: is_euler_pseudoprime(3, 0))
  23. raises(ValueError, lambda: is_euler_pseudoprime(15, 6))
  24. # A006970
  25. euler_prp = [341, 561, 1105, 1729, 1905, 2047, 2465, 3277,
  26. 4033, 4681, 5461, 6601, 8321, 8481, 10261, 10585]
  27. for p in euler_prp:
  28. assert is_euler_pseudoprime(p, 2)
  29. # A048950
  30. euler_prp = [121, 703, 1729, 1891, 2821, 3281, 7381, 8401, 8911, 10585,
  31. 12403, 15457, 15841, 16531, 18721, 19345, 23521, 24661, 28009]
  32. for p in euler_prp:
  33. assert is_euler_pseudoprime(p, 3)
  34. # A033181
  35. absolute_euler_prp = [1729, 2465, 15841, 41041, 46657, 75361,
  36. 162401, 172081, 399001, 449065, 488881]
  37. for p in absolute_euler_prp:
  38. for a in range(2, p):
  39. if gcd(a, p) != 1:
  40. continue
  41. assert is_euler_pseudoprime(p, a)
  42. def test_is_euler_jacobi_pseudoprime():
  43. assert is_euler_jacobi_pseudoprime(11, 1)
  44. assert is_euler_jacobi_pseudoprime(15, 1)
  45. def test_lucas_extrastrong_params():
  46. assert _lucas_extrastrong_params(3) == (5, 3, 1)
  47. assert _lucas_extrastrong_params(5) == (12, 4, 1)
  48. assert _lucas_extrastrong_params(7) == (5, 3, 1)
  49. assert _lucas_extrastrong_params(9) == (0, 0, 0)
  50. assert _lucas_extrastrong_params(11) == (21, 5, 1)
  51. assert _lucas_extrastrong_params(59) == (32, 6, 1)
  52. assert _lucas_extrastrong_params(479) == (117, 11, 1)
  53. def test_is_extra_strong_lucas_prp():
  54. assert is_extra_strong_lucas_prp(4) == False
  55. assert is_extra_strong_lucas_prp(989) == True
  56. assert is_extra_strong_lucas_prp(10877) == True
  57. assert is_extra_strong_lucas_prp(9) == False
  58. assert is_extra_strong_lucas_prp(16) == False
  59. assert is_extra_strong_lucas_prp(169) == False
  60. @slow
  61. def test_prps():
  62. oddcomposites = [n for n in range(1, 10**5) if
  63. n % 2 and not isprime(n)]
  64. # A checksum would be better.
  65. assert sum(oddcomposites) == 2045603465
  66. assert [n for n in oddcomposites if mr(n, [2])] == [
  67. 2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141,
  68. 52633, 65281, 74665, 80581, 85489, 88357, 90751]
  69. assert [n for n in oddcomposites if mr(n, [3])] == [
  70. 121, 703, 1891, 3281, 8401, 8911, 10585, 12403, 16531,
  71. 18721, 19345, 23521, 31621, 44287, 47197, 55969, 63139,
  72. 74593, 79003, 82513, 87913, 88573, 97567]
  73. assert [n for n in oddcomposites if mr(n, [325])] == [
  74. 9, 25, 27, 49, 65, 81, 325, 341, 343, 697, 1141, 2059,
  75. 2149, 3097, 3537, 4033, 4681, 4941, 5833, 6517, 7987, 8911,
  76. 12403, 12913, 15043, 16021, 20017, 22261, 23221, 24649,
  77. 24929, 31841, 35371, 38503, 43213, 44173, 47197, 50041,
  78. 55909, 56033, 58969, 59089, 61337, 65441, 68823, 72641,
  79. 76793, 78409, 85879]
  80. assert not any(mr(n, [9345883071009581737]) for n in oddcomposites)
  81. assert [n for n in oddcomposites if is_lucas_prp(n)] == [
  82. 323, 377, 1159, 1829, 3827, 5459, 5777, 9071, 9179, 10877,
  83. 11419, 11663, 13919, 14839, 16109, 16211, 18407, 18971,
  84. 19043, 22499, 23407, 24569, 25199, 25877, 26069, 27323,
  85. 32759, 34943, 35207, 39059, 39203, 39689, 40309, 44099,
  86. 46979, 47879, 50183, 51983, 53663, 56279, 58519, 60377,
  87. 63881, 69509, 72389, 73919, 75077, 77219, 79547, 79799,
  88. 82983, 84419, 86063, 90287, 94667, 97019, 97439]
  89. assert [n for n in oddcomposites if is_strong_lucas_prp(n)] == [
  90. 5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309,
  91. 58519, 75077, 97439]
  92. assert [n for n in oddcomposites if is_extra_strong_lucas_prp(n)
  93. ] == [
  94. 989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059,
  95. 72389, 73919, 75077]
  96. def test_proth_test():
  97. # Proth number
  98. A080075 = [3, 5, 9, 13, 17, 25, 33, 41, 49, 57, 65,
  99. 81, 97, 113, 129, 145, 161, 177, 193]
  100. # Proth prime
  101. A080076 = [3, 5, 13, 17, 41, 97, 113, 193]
  102. for n in range(200):
  103. if n in A080075:
  104. assert proth_test(n) == (n in A080076)
  105. else:
  106. raises(ValueError, lambda: proth_test(n))
  107. def test_lucas_lehmer_primality_test():
  108. for p in sieve.primerange(3, 100):
  109. assert _lucas_lehmer_primality_test(p) == (p in MERSENNE_PRIME_EXPONENTS)
  110. def test_is_mersenne_prime():
  111. assert is_mersenne_prime(-3) is False
  112. assert is_mersenne_prime(3) is True
  113. assert is_mersenne_prime(10) is False
  114. assert is_mersenne_prime(127) is True
  115. assert is_mersenne_prime(511) is False
  116. assert is_mersenne_prime(131071) is True
  117. assert is_mersenne_prime(2147483647) is True
  118. def test_isprime():
  119. s = Sieve()
  120. s.extend(100000)
  121. ps = set(s.primerange(2, 100001))
  122. for n in range(100001):
  123. # if (n in ps) != isprime(n): print n
  124. assert (n in ps) == isprime(n)
  125. assert isprime(179424673)
  126. assert isprime(20678048681)
  127. assert isprime(1968188556461)
  128. assert isprime(2614941710599)
  129. assert isprime(65635624165761929287)
  130. assert isprime(1162566711635022452267983)
  131. assert isprime(77123077103005189615466924501)
  132. assert isprime(3991617775553178702574451996736229)
  133. assert isprime(273952953553395851092382714516720001799)
  134. assert isprime(int('''
  135. 531137992816767098689588206552468627329593117727031923199444138200403\
  136. 559860852242739162502265229285668889329486246501015346579337652707239\
  137. 409519978766587351943831270835393219031728127'''))
  138. # Some Mersenne primes
  139. assert isprime(2**61 - 1)
  140. assert isprime(2**89 - 1)
  141. assert isprime(2**607 - 1)
  142. # (but not all Mersenne's are primes
  143. assert not isprime(2**601 - 1)
  144. # pseudoprimes
  145. #-------------
  146. # to some small bases
  147. assert not isprime(2152302898747)
  148. assert not isprime(3474749660383)
  149. assert not isprime(341550071728321)
  150. assert not isprime(3825123056546413051)
  151. # passes the base set [2, 3, 7, 61, 24251]
  152. assert not isprime(9188353522314541)
  153. # large examples
  154. assert not isprime(877777777777777777777777)
  155. # conjectured psi_12 given at http://mathworld.wolfram.com/StrongPseudoprime.html
  156. assert not isprime(318665857834031151167461)
  157. # conjectured psi_17 given at http://mathworld.wolfram.com/StrongPseudoprime.html
  158. assert not isprime(564132928021909221014087501701)
  159. # Arnault's 1993 number; a factor of it is
  160. # 400958216639499605418306452084546853005188166041132508774506\
  161. # 204738003217070119624271622319159721973358216316508535816696\
  162. # 9145233813917169287527980445796800452592031836601
  163. assert not isprime(int('''
  164. 803837457453639491257079614341942108138837688287558145837488917522297\
  165. 427376533365218650233616396004545791504202360320876656996676098728404\
  166. 396540823292873879185086916685732826776177102938969773947016708230428\
  167. 687109997439976544144845341155872450633409279022275296229414984230688\
  168. 1685404326457534018329786111298960644845216191652872597534901'''))
  169. # Arnault's 1995 number; can be factored as
  170. # p1*(313*(p1 - 1) + 1)*(353*(p1 - 1) + 1) where p1 is
  171. # 296744956686855105501541746429053327307719917998530433509950\
  172. # 755312768387531717701995942385964281211880336647542183455624\
  173. # 93168782883
  174. assert not isprime(int('''
  175. 288714823805077121267142959713039399197760945927972270092651602419743\
  176. 230379915273311632898314463922594197780311092934965557841894944174093\
  177. 380561511397999942154241693397290542371100275104208013496673175515285\
  178. 922696291677532547504444585610194940420003990443211677661994962953925\
  179. 045269871932907037356403227370127845389912612030924484149472897688540\
  180. 6024976768122077071687938121709811322297802059565867'''))
  181. sieve.extend(3000)
  182. assert isprime(2819)
  183. assert not isprime(2931)
  184. raises(ValueError, lambda: isprime(2.0))
  185. raises(ValueError, lambda: isprime(Float(2)))
  186. def test_is_square():
  187. assert [i for i in range(25) if is_square(i)] == [0, 1, 4, 9, 16]
  188. # issue #17044
  189. assert not is_square(60 ** 3)
  190. assert not is_square(60 ** 5)
  191. assert not is_square(84 ** 7)
  192. assert not is_square(105 ** 9)
  193. assert not is_square(120 ** 3)
  194. def test_is_gaussianprime():
  195. assert is_gaussian_prime(7*I)
  196. assert is_gaussian_prime(7)
  197. assert is_gaussian_prime(2 + 3*I)
  198. assert not is_gaussian_prime(2 + 2*I)
  199. def test_issue_27145():
  200. #https://github.com/sympy/sympy/issues/27145
  201. assert [mr(i,[2,3,5,7]) for i in (1, 2, 6)] == [False, True, False]