Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/primetest.py: 9%

1"""

2Primality testing

4"""

6from sympy.core.numbers import igcd

7from sympy.core.power import integer_nthroot

8from sympy.core.sympify import sympify

9from sympy.external.gmpy import HAS_GMPY

10from sympy.utilities.misc import as_int

12from mpmath.libmp import bitcount as _bitlength

15def _int_tuple(*i):

16 return tuple(int(_) for _ in i)

19def is_euler_pseudoprime(n, b):

20 """Returns True if n is prime or an Euler pseudoprime to base b, else False.

22 Euler Pseudoprime : In arithmetic, an odd composite integer n is called an

23 euler pseudoprime to base a, if a and n are coprime and satisfy the modular

24 arithmetic congruence relation :

26 a ^ (n-1)/2 = + 1(mod n) or

27 a ^ (n-1)/2 = - 1(mod n)

29 (where mod refers to the modulo operation).

31 Examples

32 ========

34 >>> from sympy.ntheory.primetest import is_euler_pseudoprime

35 >>> is_euler_pseudoprime(2, 5)

36 True

38 References

39 ==========

41 .. [1] https://en.wikipedia.org/wiki/Euler_pseudoprime

42 """

43 from sympy.ntheory.factor_ import trailing

45 if not mr(n, [b]):

46 return False

48 n = as_int(n)

49 r = n - 1

50 c = pow(b, r >> trailing(r), n)

52 if c == 1:

53 return True

55 while True:

56 if c == n - 1:

57 return True

58 c = pow(c, 2, n)

59 if c == 1:

60 return False

63def is_square(n, prep=True):

64 """Return True if n == a * a for some integer a, else False.

65 If n is suspected of *not* being a square then this is a

66 quick method of confirming that it is not.

68 Examples

69 ========

71 >>> from sympy.ntheory.primetest import is_square

72 >>> is_square(25)

73 True

74 >>> is_square(2)

75 False

77 References

78 ==========

80 .. [1] https://mersenneforum.org/showpost.php?p=110896

82 See Also

83 ========

84 sympy.core.power.integer_nthroot

85 """

86 if prep:

87 n = as_int(n)

88 if n < 0:

89 return False

90 if n in (0, 1):

91 return True

92 # def magic(n):

93 # s = {x**2 % n for x in range(n)}

94 # return sum(1 << bit for bit in s)

95 # >>> print(hex(magic(128)))

96 # 0x2020212020202130202021202030213

97 # >>> print(hex(magic(99)))

98 # 0x209060049048220348a410213

99 # >>> print(hex(magic(91)))

100 # 0x102e403012a0c9862c14213

101 # >>> print(hex(magic(85)))

102 # 0x121065188e001c46298213

103 if not 0x2020212020202130202021202030213 & (1 << (n & 127)):

104 return False # e.g. 2, 3

105 m = n % (99 * 91 * 85)

106 if not 0x209060049048220348a410213 & (1 << (m % 99)):

107 return False # e.g. 17, 68

108 if not 0x102e403012a0c9862c14213 & (1 << (m % 91)):

109 return False # e.g. 97, 388

110 if not 0x121065188e001c46298213 & (1 << (m % 85)):

111 return False # e.g. 793, 1408

112 # n is either:

113 # a) odd = 4*even + 1 (and square if even = k*(k + 1))

114 # b) even with

115 # odd multiplicity of 2 --> not square, e.g. 39040

116 # even multiplicity of 2, e.g. 4, 16, 36, ..., 16324

117 # removal of factors of 2 to give an odd, and rejection if

118 # any(i%2 for i in divmod(odd - 1, 4))

119 # will give an odd number in form 4*even + 1.

120 # Use of `trailing` to check the power of 2 is not done since it

121 # does not apply to a large percentage of arbitrary numbers

122 # and the integer_nthroot is able to quickly resolve these cases.

123 return integer_nthroot(n, 2)[1]

124

125

126def _test(n, base, s, t):

127 """Miller-Rabin strong pseudoprime test for one base.

128 Return False if n is definitely composite, True if n is

129 probably prime, with a probability greater than 3/4.

130

131 """

132 # do the Fermat test

133 b = pow(base, t, n)

134 if b == 1 or b == n - 1:

135 return True

136 else:

137 for j in range(1, s):

138 b = pow(b, 2, n)

139 if b == n - 1:

140 return True

141 # see I. Niven et al. "An Introduction to Theory of Numbers", page 78

142 if b == 1:

143 return False

144 return False

145

146

147def mr(n, bases):

148 """Perform a Miller-Rabin strong pseudoprime test on n using a

149 given list of bases/witnesses.

150

151 References

152 ==========

153

154 .. [1] Richard Crandall & Carl Pomerance (2005), "Prime Numbers:

155 A Computational Perspective", Springer, 2nd edition, 135-138

156

157 A list of thresholds and the bases they require are here:

158 https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants

159

160 Examples

161 ========

162

163 >>> from sympy.ntheory.primetest import mr

164 >>> mr(1373651, [2, 3])

165 False

166 >>> mr(479001599, [31, 73])

167 True

168

169 """

170 from sympy.ntheory.factor_ import trailing

171 from sympy.polys.domains import ZZ

172

173 n = as_int(n)

174 if n < 2:

175 return False

176 # remove powers of 2 from n-1 (= t * 2**s)

177 s = trailing(n - 1)

178 t = n >> s

179 for base in bases:

180 # Bases >= n are wrapped, bases < 2 are invalid

181 if base >= n:

182 base %= n

183 if base >= 2:

184 base = ZZ(base)

185 if not _test(n, base, s, t):

186 return False

187 return True

188

189

190def _lucas_sequence(n, P, Q, k):

191 """Return the modular Lucas sequence (U_k, V_k, Q_k).

192

193 Given a Lucas sequence defined by P, Q, returns the kth values for

194 U and V, along with Q^k, all modulo n. This is intended for use with

195 possibly very large values of n and k, where the combinatorial functions

196 would be completely unusable.

197

198 The modular Lucas sequences are used in numerous places in number theory,

199 especially in the Lucas compositeness tests and the various n + 1 proofs.

200

201 Examples

202 ========

203

204 >>> from sympy.ntheory.primetest import _lucas_sequence

205 >>> N = 10**2000 + 4561

206 >>> sol = U, V, Qk = _lucas_sequence(N, 3, 1, N//2); sol

207 (0, 2, 1)

208

209 """

210 D = P*P - 4*Q

211 if n < 2:

212 raise ValueError("n must be >= 2")

213 if k < 0:

214 raise ValueError("k must be >= 0")

215 if D == 0:

216 raise ValueError("D must not be zero")

217

218 if k == 0:

219 return _int_tuple(0, 2, Q)

220 U = 1

221 V = P

222 Qk = Q

223 b = _bitlength(k)

224 if Q == 1:

225 # Optimization for extra strong tests.

226 while b > 1:

227 U = (U*V) % n

228 V = (V*V - 2) % n

229 b -= 1

230 if (k >> (b - 1)) & 1:

231 U, V = U*P + V, V*P + U*D

232 if U & 1:

233 U += n

234 if V & 1:

235 V += n

236 U, V = U >> 1, V >> 1

237 elif P == 1 and Q == -1:

238 # Small optimization for 50% of Selfridge parameters.

239 while b > 1:

240 U = (U*V) % n

241 if Qk == 1:

242 V = (V*V - 2) % n

243 else:

244 V = (V*V + 2) % n

245 Qk = 1

246 b -= 1

247 if (k >> (b-1)) & 1:

248 U, V = U + V, V + U*D

249 if U & 1:

250 U += n

251 if V & 1:

252 V += n

253 U, V = U >> 1, V >> 1

254 Qk = -1

255 else:

256 # The general case with any P and Q.

257 while b > 1:

258 U = (U*V) % n

259 V = (V*V - 2*Qk) % n

260 Qk *= Qk

261 b -= 1

262 if (k >> (b - 1)) & 1:

263 U, V = U*P + V, V*P + U*D

264 if U & 1:

265 U += n

266 if V & 1:

267 V += n

268 U, V = U >> 1, V >> 1

269 Qk *= Q

270 Qk %= n

271 return _int_tuple(U % n, V % n, Qk)

272

273

274def _lucas_selfridge_params(n):

275 """Calculates the Selfridge parameters (D, P, Q) for n. This is

276 method A from page 1401 of Baillie and Wagstaff.

277

278 References

279 ==========

280 .. [1] "Lucas Pseudoprimes", Baillie and Wagstaff, 1980.

281 http://mpqs.free.fr/LucasPseudoprimes.pdf

282 """

283 from sympy.ntheory.residue_ntheory import jacobi_symbol

284 D = 5

285 while True:

286 g = igcd(abs(D), n)

287 if g > 1 and g != n:

288 return (0, 0, 0)

289 if jacobi_symbol(D, n) == -1:

290 break

291 if D > 0:

292 D = -D - 2

293 else:

294 D = -D + 2

295 return _int_tuple(D, 1, (1 - D)/4)

296

297

298def _lucas_extrastrong_params(n):

299 """Calculates the "extra strong" parameters (D, P, Q) for n.

300

301 References

302 ==========

303 .. [1] OEIS A217719: Extra Strong Lucas Pseudoprimes

304 https://oeis.org/A217719

305 .. [1] https://en.wikipedia.org/wiki/Lucas_pseudoprime

306 """

307 from sympy.ntheory.residue_ntheory import jacobi_symbol

308 P, Q, D = 3, 1, 5

309 while True:

310 g = igcd(D, n)

311 if g > 1 and g != n:

312 return (0, 0, 0)

313 if jacobi_symbol(D, n) == -1:

314 break

315 P += 1

316 D = P*P - 4

317 return _int_tuple(D, P, Q)

318

319

320def is_lucas_prp(n):

321 """Standard Lucas compositeness test with Selfridge parameters. Returns

322 False if n is definitely composite, and True if n is a Lucas probable

323 prime.

324

325 This is typically used in combination with the Miller-Rabin test.

326

327 References

328 ==========

329 - "Lucas Pseudoprimes", Baillie and Wagstaff, 1980.

330 http://mpqs.free.fr/LucasPseudoprimes.pdf

331 - OEIS A217120: Lucas Pseudoprimes

332 https://oeis.org/A217120

333 - https://en.wikipedia.org/wiki/Lucas_pseudoprime

334

335 Examples

336 ========

337

338 >>> from sympy.ntheory.primetest import isprime, is_lucas_prp

339 >>> for i in range(10000):

340 ... if is_lucas_prp(i) and not isprime(i):

341 ... print(i)

342 323

343 377

344 1159

345 1829

346 3827

347 5459

348 5777

349 9071

350 9179

351 """

352 n = as_int(n)

353 if n == 2:

354 return True

355 if n < 2 or (n % 2) == 0:

356 return False

357 if is_square(n, False):

358 return False

359

360 D, P, Q = _lucas_selfridge_params(n)

361 if D == 0:

362 return False

363 U, V, Qk = _lucas_sequence(n, P, Q, n+1)

364 return U == 0

365

366

367def is_strong_lucas_prp(n):

368 """Strong Lucas compositeness test with Selfridge parameters. Returns

369 False if n is definitely composite, and True if n is a strong Lucas

370 probable prime.

371

372 This is often used in combination with the Miller-Rabin test, and

373 in particular, when combined with M-R base 2 creates the strong BPSW test.

374

375 References

376 ==========

377 - "Lucas Pseudoprimes", Baillie and Wagstaff, 1980.

378 http://mpqs.free.fr/LucasPseudoprimes.pdf

379 - OEIS A217255: Strong Lucas Pseudoprimes

380 https://oeis.org/A217255

381 - https://en.wikipedia.org/wiki/Lucas_pseudoprime

382 - https://en.wikipedia.org/wiki/Baillie-PSW_primality_test

383

384 Examples

385 ========

386

387 >>> from sympy.ntheory.primetest import isprime, is_strong_lucas_prp

388 >>> for i in range(20000):

389 ... if is_strong_lucas_prp(i) and not isprime(i):

390 ... print(i)

391 5459

392 5777

393 10877

394 16109

395 18971

396 """

397 from sympy.ntheory.factor_ import trailing

398 n = as_int(n)

399 if n == 2:

400 return True

401 if n < 2 or (n % 2) == 0:

402 return False

403 if is_square(n, False):

404 return False

405

406 D, P, Q = _lucas_selfridge_params(n)

407 if D == 0:

408 return False

409

410 # remove powers of 2 from n+1 (= k * 2**s)

411 s = trailing(n + 1)

412 k = (n+1) >> s

413

414 U, V, Qk = _lucas_sequence(n, P, Q, k)

415

416 if U == 0 or V == 0:

417 return True

418 for r in range(1, s):

419 V = (V*V - 2*Qk) % n

420 if V == 0:

421 return True

422 Qk = pow(Qk, 2, n)

423 return False

424

425

426def is_extra_strong_lucas_prp(n):

427 """Extra Strong Lucas compositeness test. Returns False if n is

428 definitely composite, and True if n is a "extra strong" Lucas probable

429 prime.

430

431 The parameters are selected using P = 3, Q = 1, then incrementing P until

432 (D|n) == -1. The test itself is as defined in Grantham 2000, from the

433 Mo and Jones preprint. The parameter selection and test are the same as

434 used in OEIS A217719, Perl's Math::Prime::Util, and the Lucas pseudoprime

435 page on Wikipedia.

436

437 With these parameters, there are no counterexamples below 2^64 nor any

438 known above that range. It is 20-50% faster than the strong test.

439

440 Because of the different parameters selected, there is no relationship

441 between the strong Lucas pseudoprimes and extra strong Lucas pseudoprimes.

442 In particular, one is not a subset of the other.

443

444 References

445 ==========

446 - "Frobenius Pseudoprimes", Jon Grantham, 2000.

447 https://www.ams.org/journals/mcom/2001-70-234/S0025-5718-00-01197-2/

448 - OEIS A217719: Extra Strong Lucas Pseudoprimes

449 https://oeis.org/A217719

450 - https://en.wikipedia.org/wiki/Lucas_pseudoprime

451

452 Examples

453 ========

454

455 >>> from sympy.ntheory.primetest import isprime, is_extra_strong_lucas_prp

456 >>> for i in range(20000):

457 ... if is_extra_strong_lucas_prp(i) and not isprime(i):

458 ... print(i)

459 989

460 3239

461 5777

462 10877

463 """

464 # Implementation notes:

465 # 1) the parameters differ from Thomas R. Nicely's. His parameter

466 # selection leads to pseudoprimes that overlap M-R tests, and

467 # contradict Baillie and Wagstaff's suggestion of (D|n) = -1.

468 # 2) The MathWorld page as of June 2013 specifies Q=-1. The Lucas

469 # sequence must have Q=1. See Grantham theorem 2.3, any of the

470 # references on the MathWorld page, or run it and see Q=-1 is wrong.

471 from sympy.ntheory.factor_ import trailing

472 n = as_int(n)

473 if n == 2:

474 return True

475 if n < 2 or (n % 2) == 0:

476 return False

477 if is_square(n, False):

478 return False

479

480 D, P, Q = _lucas_extrastrong_params(n)

481 if D == 0:

482 return False

483

484 # remove powers of 2 from n+1 (= k * 2**s)

485 s = trailing(n + 1)

486 k = (n+1) >> s

487

488 U, V, Qk = _lucas_sequence(n, P, Q, k)

489

490 if U == 0 and (V == 2 or V == n - 2):

491 return True

492 for r in range(1, s):

493 if V == 0:

494 return True

495 V = (V*V - 2) % n

496 return False

497

498

499def isprime(n):

500 """

501 Test if n is a prime number (True) or not (False). For n < 2^64 the

502 answer is definitive; larger n values have a small probability of actually

503 being pseudoprimes.

504

505 Negative numbers (e.g. -2) are not considered prime.

506

507 The first step is looking for trivial factors, which if found enables

508 a quick return. Next, if the sieve is large enough, use bisection search

509 on the sieve. For small numbers, a set of deterministic Miller-Rabin

510 tests are performed with bases that are known to have no counterexamples

511 in their range. Finally if the number is larger than 2^64, a strong

512 BPSW test is performed. While this is a probable prime test and we

513 believe counterexamples exist, there are no known counterexamples.

514

515 Examples

516 ========

517

518 >>> from sympy.ntheory import isprime

519 >>> isprime(13)

520 True

521 >>> isprime(13.0) # limited precision

522 False

523 >>> isprime(15)

524 False

525

526 Notes

527 =====

528

529 This routine is intended only for integer input, not numerical

530 expressions which may represent numbers. Floats are also

531 rejected as input because they represent numbers of limited

532 precision. While it is tempting to permit 7.0 to represent an

533 integer there are errors that may "pass silently" if this is

534 allowed:

535

536 >>> from sympy import Float, S

537 >>> int(1e3) == 1e3 == 10**3

538 True

539 >>> int(1e23) == 1e23

540 True

541 >>> int(1e23) == 10**23

542 False

543

544 >>> near_int = 1 + S(1)/10**19

545 >>> near_int == int(near_int)

546 False

547 >>> n = Float(near_int, 10) # truncated by precision

548 >>> n == int(n)

549 True

550 >>> n = Float(near_int, 20)

551 >>> n == int(n)

552 False

553

554 See Also

555 ========

556

557 sympy.ntheory.generate.primerange : Generates all primes in a given range

558 sympy.ntheory.generate.primepi : Return the number of primes less than or equal to n

559 sympy.ntheory.generate.prime : Return the nth prime

560

561 References

562 ==========

563 - https://en.wikipedia.org/wiki/Strong_pseudoprime

564 - "Lucas Pseudoprimes", Baillie and Wagstaff, 1980.

565 http://mpqs.free.fr/LucasPseudoprimes.pdf

566 - https://en.wikipedia.org/wiki/Baillie-PSW_primality_test

567 """

568 try:

569 n = as_int(n)

570 except ValueError:

571 return False

572

573 # Step 1, do quick composite testing via trial division. The individual

574 # modulo tests benchmark faster than one or two primorial igcds for me.

575 # The point here is just to speedily handle small numbers and many

576 # composites. Step 2 only requires that n <= 2 get handled here.

577 if n in [2, 3, 5]:

578 return True

579 if n < 2 or (n % 2) == 0 or (n % 3) == 0 or (n % 5) == 0:

580 return False

581 if n < 49:

582 return True

583 if (n % 7) == 0 or (n % 11) == 0 or (n % 13) == 0 or (n % 17) == 0 or \

584 (n % 19) == 0 or (n % 23) == 0 or (n % 29) == 0 or (n % 31) == 0 or \

585 (n % 37) == 0 or (n % 41) == 0 or (n % 43) == 0 or (n % 47) == 0:

586 return False

587 if n < 2809:

588 return True

589 if n < 31417:

590 return pow(2, n, n) == 2 and n not in [7957, 8321, 13747, 18721, 19951, 23377]

591

592 # bisection search on the sieve if the sieve is large enough

593 from sympy.ntheory.generate import sieve as s

594 if n <= s._list[-1]:

595 l, u = s.search(n)

596 return l == u

597

598 # If we have GMPY2, skip straight to step 3 and do a strong BPSW test.

599 # This should be a bit faster than our step 2, and for large values will

600 # be a lot faster than our step 3 (C+GMP vs. Python).

601 if HAS_GMPY == 2:

602 from gmpy2 import is_strong_prp, is_strong_selfridge_prp

603 return is_strong_prp(n, 2) and is_strong_selfridge_prp(n)

604

605

606 # Step 2: deterministic Miller-Rabin testing for numbers < 2^64. See:

607 # https://miller-rabin.appspot.com/

608 # for lists. We have made sure the M-R routine will successfully handle

609 # bases larger than n, so we can use the minimal set.

610 # In September 2015 deterministic numbers were extended to over 2^81.

611 # https://arxiv.org/pdf/1509.00864.pdf

612 # https://oeis.org/A014233

613 if n < 341531:

614 return mr(n, [9345883071009581737])

615 if n < 885594169:

616 return mr(n, [725270293939359937, 3569819667048198375])

617 if n < 350269456337:

618 return mr(n, [4230279247111683200, 14694767155120705706, 16641139526367750375])

619 if n < 55245642489451:

620 return mr(n, [2, 141889084524735, 1199124725622454117, 11096072698276303650])

621 if n < 7999252175582851:

622 return mr(n, [2, 4130806001517, 149795463772692060, 186635894390467037, 3967304179347715805])

623 if n < 585226005592931977:

624 return mr(n, [2, 123635709730000, 9233062284813009, 43835965440333360, 761179012939631437, 1263739024124850375])

625 if n < 18446744073709551616:

626 return mr(n, [2, 325, 9375, 28178, 450775, 9780504, 1795265022])

627 if n < 318665857834031151167461:

628 return mr(n, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37])

629 if n < 3317044064679887385961981:

630 return mr(n, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41])

631

632 # We could do this instead at any point:

633 #if n < 18446744073709551616:

634 # return mr(n, [2]) and is_extra_strong_lucas_prp(n)

635

636 # Here are tests that are safe for MR routines that don't understand

637 # large bases.

638 #if n < 9080191:

639 # return mr(n, [31, 73])

640 #if n < 19471033:

641 # return mr(n, [2, 299417])

642 #if n < 38010307:

643 # return mr(n, [2, 9332593])

644 #if n < 316349281:

645 # return mr(n, [11000544, 31481107])

646 #if n < 4759123141:

647 # return mr(n, [2, 7, 61])

648 #if n < 105936894253:

649 # return mr(n, [2, 1005905886, 1340600841])

650 #if n < 31858317218647:

651 # return mr(n, [2, 642735, 553174392, 3046413974])

652 #if n < 3071837692357849:

653 # return mr(n, [2, 75088, 642735, 203659041, 3613982119])

654 #if n < 18446744073709551616:

655 # return mr(n, [2, 325, 9375, 28178, 450775, 9780504, 1795265022])

656

657 # Step 3: BPSW.

658 #

659 # Time for isprime(10**2000 + 4561), no gmpy or gmpy2 installed

660 # 44.0s old isprime using 46 bases

661 # 5.3s strong BPSW + one random base

662 # 4.3s extra strong BPSW + one random base

663 # 4.1s strong BPSW

664 # 3.2s extra strong BPSW

665

666 # Classic BPSW from page 1401 of the paper. See alternate ideas below.

667 return mr(n, [2]) and is_strong_lucas_prp(n)

668

669 # Using extra strong test, which is somewhat faster

670 #return mr(n, [2]) and is_extra_strong_lucas_prp(n)

671

672 # Add a random M-R base

673 #import random

674 #return mr(n, [2, random.randint(3, n-1)]) and is_strong_lucas_prp(n)

675

676

677def is_gaussian_prime(num):

678 r"""Test if num is a Gaussian prime number.

679

680 References

681 ==========

682

683 .. [1] https://oeis.org/wiki/Gaussian_primes

684 """

685

686 num = sympify(num)

687 a, b = num.as_real_imag()

688 a = as_int(a, strict=False)

689 b = as_int(b, strict=False)

690 if a == 0:

691 b = abs(b)

692 return isprime(b) and b % 4 == 3

693 elif b == 0:

694 a = abs(a)

695 return isprime(a) and a % 4 == 3

696 return isprime(a**2 + b**2)