Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/generate.py: 16%

1""" 

2Generating and counting primes. 

3 

4""" 

5 

6import random 

7from bisect import bisect 

8from itertools import count 

9# Using arrays for sieving instead of lists greatly reduces 

10# memory consumption 

11from array import array as _array 

12 

13from sympy.core.function import Function 

14from sympy.core.singleton import S 

15from .primetest import isprime 

16from sympy.utilities.misc import as_int 

17 

18 

19def _azeros(n): 

20 return _array('l', [0]*n) 

21 

22 

23def _aset(*v): 

24 return _array('l', v) 

25 

26 

27def _arange(a, b): 

28 return _array('l', range(a, b)) 

29 

30 

31def _as_int_ceiling(a): 

32 """ Wrapping ceiling in as_int will raise an error if there was a problem 

33 determining whether the expression was exactly an integer or not.""" 

34 from sympy.functions.elementary.integers import ceiling 

35 return as_int(ceiling(a)) 

36 

37 

38class Sieve: 

39 """An infinite list of prime numbers, implemented as a dynamically 

40 growing sieve of Eratosthenes. When a lookup is requested involving 

41 an odd number that has not been sieved, the sieve is automatically 

42 extended up to that number. 

43 

44 Examples 

45 ======== 

46 

47 >>> from sympy import sieve 

48 >>> sieve._reset() # this line for doctest only 

49 >>> 25 in sieve 

50 False 

51 >>> sieve._list 

52 array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23]) 

53 """ 

54 

55 # data shared (and updated) by all Sieve instances 

56 def __init__(self): 

57 self._n = 6 

58 self._list = _aset(2, 3, 5, 7, 11, 13) # primes 

59 self._tlist = _aset(0, 1, 1, 2, 2, 4) # totient 

60 self._mlist = _aset(0, 1, -1, -1, 0, -1) # mobius 

61 assert all(len(i) == self._n for i in (self._list, self._tlist, self._mlist)) 

62 

63 def __repr__(self): 

64 return ("<%s sieve (%i): %i, %i, %i, ... %i, %i\n" 

65 "%s sieve (%i): %i, %i, %i, ... %i, %i\n" 

66 "%s sieve (%i): %i, %i, %i, ... %i, %i>") % ( 

67 'prime', len(self._list), 

68 self._list[0], self._list[1], self._list[2], 

69 self._list[-2], self._list[-1], 

70 'totient', len(self._tlist), 

71 self._tlist[0], self._tlist[1], 

72 self._tlist[2], self._tlist[-2], self._tlist[-1], 

73 'mobius', len(self._mlist), 

74 self._mlist[0], self._mlist[1], 

75 self._mlist[2], self._mlist[-2], self._mlist[-1]) 

76 

77 def _reset(self, prime=None, totient=None, mobius=None): 

78 """Reset all caches (default). To reset one or more set the 

79 desired keyword to True.""" 

80 if all(i is None for i in (prime, totient, mobius)): 

81 prime = totient = mobius = True 

82 if prime: 

83 self._list = self._list[:self._n] 

84 if totient: 

85 self._tlist = self._tlist[:self._n] 

86 if mobius: 

87 self._mlist = self._mlist[:self._n] 

88 

89 def extend(self, n): 

90 """Grow the sieve to cover all primes <= n (a real number). 

91 

92 Examples 

93 ======== 

94 

95 >>> from sympy import sieve 

96 >>> sieve._reset() # this line for doctest only 

97 >>> sieve.extend(30) 

98 >>> sieve[10] == 29 

99 True 

100 """ 

101 n = int(n) 

102 if n <= self._list[-1]: 

103 return 

104 

105 # We need to sieve against all bases up to sqrt(n). 

106 # This is a recursive call that will do nothing if there are enough 

107 # known bases already. 

108 maxbase = int(n**0.5) + 1 

109 self.extend(maxbase) 

110 

111 # Create a new sieve starting from sqrt(n) 

112 begin = self._list[-1] + 1 

113 newsieve = _arange(begin, n + 1) 

114 

115 # Now eliminate all multiples of primes in [2, sqrt(n)] 

116 for p in self.primerange(maxbase): 

117 # Start counting at a multiple of p, offsetting 

118 # the index to account for the new sieve's base index 

119 startindex = (-begin) % p 

120 for i in range(startindex, len(newsieve), p): 

121 newsieve[i] = 0 

122 

123 # Merge the sieves 

124 self._list += _array('l', [x for x in newsieve if x]) 

125 

126 def extend_to_no(self, i): 

127 """Extend to include the ith prime number. 

128 

129 Parameters 

130 ========== 

131 

132 i : integer 

133 

134 Examples 

135 ======== 

136 

137 >>> from sympy import sieve 

138 >>> sieve._reset() # this line for doctest only 

139 >>> sieve.extend_to_no(9) 

140 >>> sieve._list 

141 array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23]) 

142 

143 Notes 

144 ===== 

145 

146 The list is extended by 50% if it is too short, so it is 

147 likely that it will be longer than requested. 

148 """ 

149 i = as_int(i) 

150 while len(self._list) < i: 

151 self.extend(int(self._list[-1] * 1.5)) 

152 

153 def primerange(self, a, b=None): 

154 """Generate all prime numbers in the range [2, a) or [a, b). 

155 

156 Examples 

157 ======== 

158 

159 >>> from sympy import sieve, prime 

160 

161 All primes less than 19: 

162 

163 >>> print([i for i in sieve.primerange(19)]) 

164 [2, 3, 5, 7, 11, 13, 17] 

165 

166 All primes greater than or equal to 7 and less than 19: 

167 

168 >>> print([i for i in sieve.primerange(7, 19)]) 

169 [7, 11, 13, 17] 

170 

171 All primes through the 10th prime 

172 

173 >>> list(sieve.primerange(prime(10) + 1)) 

174 [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] 

175 

176 """ 

177 

178 if b is None: 

179 b = _as_int_ceiling(a) 

180 a = 2 

181 else: 

182 a = max(2, _as_int_ceiling(a)) 

183 b = _as_int_ceiling(b) 

184 if a >= b: 

185 return 

186 self.extend(b) 

187 i = self.search(a)[1] 

188 maxi = len(self._list) + 1 

189 while i < maxi: 

190 p = self._list[i - 1] 

191 if p < b: 

192 yield p 

193 i += 1 

194 else: 

195 return 

196 

197 def totientrange(self, a, b): 

198 """Generate all totient numbers for the range [a, b). 

199 

200 Examples 

201 ======== 

202 

203 >>> from sympy import sieve 

204 >>> print([i for i in sieve.totientrange(7, 18)]) 

205 [6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16] 

206 """ 

207 a = max(1, _as_int_ceiling(a)) 

208 b = _as_int_ceiling(b) 

209 n = len(self._tlist) 

210 if a >= b: 

211 return 

212 elif b <= n: 

213 for i in range(a, b): 

214 yield self._tlist[i] 

215 else: 

216 self._tlist += _arange(n, b) 

217 for i in range(1, n): 

218 ti = self._tlist[i] 

219 startindex = (n + i - 1) // i * i 

220 for j in range(startindex, b, i): 

221 self._tlist[j] -= ti 

222 if i >= a: 

223 yield ti 

224 

225 for i in range(n, b): 

226 ti = self._tlist[i] 

227 for j in range(2 * i, b, i): 

228 self._tlist[j] -= ti 

229 if i >= a: 

230 yield ti 

231 

232 def mobiusrange(self, a, b): 

233 """Generate all mobius numbers for the range [a, b). 

234 

235 Parameters 

236 ========== 

237 

238 a : integer 

239 First number in range 

240 

241 b : integer 

242 First number outside of range 

243 

244 Examples 

245 ======== 

246 

247 >>> from sympy import sieve 

248 >>> print([i for i in sieve.mobiusrange(7, 18)]) 

249 [-1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1] 

250 """ 

251 a = max(1, _as_int_ceiling(a)) 

252 b = _as_int_ceiling(b) 

253 n = len(self._mlist) 

254 if a >= b: 

255 return 

256 elif b <= n: 

257 for i in range(a, b): 

258 yield self._mlist[i] 

259 else: 

260 self._mlist += _azeros(b - n) 

261 for i in range(1, n): 

262 mi = self._mlist[i] 

263 startindex = (n + i - 1) // i * i 

264 for j in range(startindex, b, i): 

265 self._mlist[j] -= mi 

266 if i >= a: 

267 yield mi 

268 

269 for i in range(n, b): 

270 mi = self._mlist[i] 

271 for j in range(2 * i, b, i): 

272 self._mlist[j] -= mi 

273 if i >= a: 

274 yield mi 

275 

276 def search(self, n): 

277 """Return the indices i, j of the primes that bound n. 

278 

279 If n is prime then i == j. 

280 

281 Although n can be an expression, if ceiling cannot convert 

282 it to an integer then an n error will be raised. 

283 

284 Examples 

285 ======== 

286 

287 >>> from sympy import sieve 

288 >>> sieve.search(25) 

289 (9, 10) 

290 >>> sieve.search(23) 

291 (9, 9) 

292 """ 

293 test = _as_int_ceiling(n) 

294 n = as_int(n) 

295 if n < 2: 

296 raise ValueError("n should be >= 2 but got: %s" % n) 

297 if n > self._list[-1]: 

298 self.extend(n) 

299 b = bisect(self._list, n) 

300 if self._list[b - 1] == test: 

301 return b, b 

302 else: 

303 return b, b + 1 

304 

305 def __contains__(self, n): 

306 try: 

307 n = as_int(n) 

308 assert n >= 2 

309 except (ValueError, AssertionError): 

310 return False 

311 if n % 2 == 0: 

312 return n == 2 

313 a, b = self.search(n) 

314 return a == b 

315 

316 def __iter__(self): 

317 for n in count(1): 

318 yield self[n] 

319 

320 def __getitem__(self, n): 

321 """Return the nth prime number""" 

322 if isinstance(n, slice): 

323 self.extend_to_no(n.stop) 

324 # Python 2.7 slices have 0 instead of None for start, so 

325 # we can't default to 1. 

326 start = n.start if n.start is not None else 0 

327 if start < 1: 

328 # sieve[:5] would be empty (starting at -1), let's 

329 # just be explicit and raise. 

330 raise IndexError("Sieve indices start at 1.") 

331 return self._list[start - 1:n.stop - 1:n.step] 

332 else: 

333 if n < 1: 

334 # offset is one, so forbid explicit access to sieve[0] 

335 # (would surprisingly return the last one). 

336 raise IndexError("Sieve indices start at 1.") 

337 n = as_int(n) 

338 self.extend_to_no(n) 

339 return self._list[n - 1] 

340 

341# Generate a global object for repeated use in trial division etc 

342sieve = Sieve() 

343 

344 

345def prime(nth): 

346 r""" Return the nth prime, with the primes indexed as prime(1) = 2, 

347 prime(2) = 3, etc.... The nth prime is approximately $n\log(n)$. 

348 

349 Logarithmic integral of $x$ is a pretty nice approximation for number of 

350 primes $\le x$, i.e. 

351 li(x) ~ pi(x) 

352 In fact, for the numbers we are concerned about( x<1e11 ), 

353 li(x) - pi(x) < 50000 

354 

355 Also, 

356 li(x) > pi(x) can be safely assumed for the numbers which 

357 can be evaluated by this function. 

358 

359 Here, we find the least integer m such that li(m) > n using binary search. 

360 Now pi(m-1) < li(m-1) <= n, 

361 

362 We find pi(m - 1) using primepi function. 

363 

364 Starting from m, we have to find n - pi(m-1) more primes. 

365 

366 For the inputs this implementation can handle, we will have to test 

367 primality for at max about 10**5 numbers, to get our answer. 

368 

369 Examples 

370 ======== 

371 

372 >>> from sympy import prime 

373 >>> prime(10) 

374 29 

375 >>> prime(1) 

376 2 

377 >>> prime(100000) 

378 1299709 

379 

380 See Also 

381 ======== 

382 

383 sympy.ntheory.primetest.isprime : Test if n is prime 

384 primerange : Generate all primes in a given range 

385 primepi : Return the number of primes less than or equal to n 

386 

387 References 

388 ========== 

389 

390 .. [1] https://en.wikipedia.org/wiki/Prime_number_theorem#Table_of_.CF.80.28x.29.2C_x_.2F_log_x.2C_and_li.28x.29 

391 .. [2] https://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number 

392 .. [3] https://en.wikipedia.org/wiki/Skewes%27_number 

393 """ 

394 n = as_int(nth) 

395 if n < 1: 

396 raise ValueError("nth must be a positive integer; prime(1) == 2") 

397 if n <= len(sieve._list): 

398 return sieve[n] 

399 

400 from sympy.functions.elementary.exponential import log 

401 from sympy.functions.special.error_functions import li 

402 a = 2 # Lower bound for binary search 

403 b = int(n*(log(n) + log(log(n)))) # Upper bound for the search. 

404 

405 while a < b: 

406 mid = (a + b) >> 1 

407 if li(mid) > n: 

408 b = mid 

409 else: 

410 a = mid + 1 

411 n_primes = primepi(a - 1) 

412 while n_primes < n: 

413 if isprime(a): 

414 n_primes += 1 

415 a += 1 

416 return a - 1 

417 

418 

419class primepi(Function): 

420 r""" Represents the prime counting function pi(n) = the number 

421 of prime numbers less than or equal to n. 

422 

423 Algorithm Description: 

424 

425 In sieve method, we remove all multiples of prime p 

426 except p itself. 

427 

428 Let phi(i,j) be the number of integers 2 <= k <= i 

429 which remain after sieving from primes less than 

430 or equal to j. 

431 Clearly, pi(n) = phi(n, sqrt(n)) 

432 

433 If j is not a prime, 

434 phi(i,j) = phi(i, j - 1) 

435 

436 if j is a prime, 

437 We remove all numbers(except j) whose 

438 smallest prime factor is j. 

439 

440 Let $x= j \times a$ be such a number, where $2 \le a \le i / j$ 

441 Now, after sieving from primes $\le j - 1$, 

442 a must remain 

443 (because x, and hence a has no prime factor $\le j - 1$) 

444 Clearly, there are phi(i / j, j - 1) such a 

445 which remain on sieving from primes $\le j - 1$ 

446 

447 Now, if a is a prime less than equal to j - 1, 

448 $x= j \times a$ has smallest prime factor = a, and 

449 has already been removed(by sieving from a). 

450 So, we do not need to remove it again. 

451 (Note: there will be pi(j - 1) such x) 

452 

453 Thus, number of x, that will be removed are: 

454 phi(i / j, j - 1) - phi(j - 1, j - 1) 

455 (Note that pi(j - 1) = phi(j - 1, j - 1)) 

456 

457 $\Rightarrow$ phi(i,j) = phi(i, j - 1) - phi(i / j, j - 1) + phi(j - 1, j - 1) 

458 

459 So,following recursion is used and implemented as dp: 

460 

461 phi(a, b) = phi(a, b - 1), if b is not a prime 

462 phi(a, b) = phi(a, b-1)-phi(a / b, b-1) + phi(b-1, b-1), if b is prime 

463 

464 Clearly a is always of the form floor(n / k), 

465 which can take at most $2\sqrt{n}$ values. 

466 Two arrays arr1,arr2 are maintained 

467 arr1[i] = phi(i, j), 

468 arr2[i] = phi(n // i, j) 

469 

470 Finally the answer is arr2[1] 

471 

472 Examples 

473 ======== 

474 

475 >>> from sympy import primepi, prime, prevprime, isprime 

476 >>> primepi(25) 

477 9 

478 

479 So there are 9 primes less than or equal to 25. Is 25 prime? 

480 

481 >>> isprime(25) 

482 False 

483 

484 It is not. So the first prime less than 25 must be the 

485 9th prime: 

486 

487 >>> prevprime(25) == prime(9) 

488 True 

489 

490 See Also 

491 ======== 

492 

493 sympy.ntheory.primetest.isprime : Test if n is prime 

494 primerange : Generate all primes in a given range 

495 prime : Return the nth prime 

496 """ 

497 @classmethod 

498 def eval(cls, n): 

499 if n is S.Infinity: 

500 return S.Infinity 

501 if n is S.NegativeInfinity: 

502 return S.Zero 

503 

504 try: 

505 n = int(n) 

506 except TypeError: 

507 if n.is_real == False or n is S.NaN: 

508 raise ValueError("n must be real") 

509 return 

510 

511 if n < 2: 

512 return S.Zero 

513 if n <= sieve._list[-1]: 

514 return S(sieve.search(n)[0]) 

515 lim = int(n ** 0.5) 

516 lim -= 1 

517 lim = max(lim, 0) 

518 while lim * lim <= n: 

519 lim += 1 

520 lim -= 1 

521 arr1 = [0] * (lim + 1) 

522 arr2 = [0] * (lim + 1) 

523 for i in range(1, lim + 1): 

524 arr1[i] = i - 1 

525 arr2[i] = n // i - 1 

526 for i in range(2, lim + 1): 

527 # Presently, arr1[k]=phi(k,i - 1), 

528 # arr2[k] = phi(n // k,i - 1) 

529 if arr1[i] == arr1[i - 1]: 

530 continue 

531 p = arr1[i - 1] 

532 for j in range(1, min(n // (i * i), lim) + 1): 

533 st = i * j 

534 if st <= lim: 

535 arr2[j] -= arr2[st] - p 

536 else: 

537 arr2[j] -= arr1[n // st] - p 

538 lim2 = min(lim, i * i - 1) 

539 for j in range(lim, lim2, -1): 

540 arr1[j] -= arr1[j // i] - p 

541 return S(arr2[1]) 

542 

543 

544def nextprime(n, ith=1): 

545 """ Return the ith prime greater than n. 

546 

547 i must be an integer. 

548 

549 Notes 

550 ===== 

551 

552 Potential primes are located at 6*j +/- 1. This 

553 property is used during searching. 

554 

555 >>> from sympy import nextprime 

556 >>> [(i, nextprime(i)) for i in range(10, 15)] 

557 [(10, 11), (11, 13), (12, 13), (13, 17), (14, 17)] 

558 >>> nextprime(2, ith=2) # the 2nd prime after 2 

559 5 

560 

561 See Also 

562 ======== 

563 

564 prevprime : Return the largest prime smaller than n 

565 primerange : Generate all primes in a given range 

566 

567 """ 

568 n = int(n) 

569 i = as_int(ith) 

570 if i > 1: 

571 pr = n 

572 j = 1 

573 while 1: 

574 pr = nextprime(pr) 

575 j += 1 

576 if j > i: 

577 break 

578 return pr 

579 

580 if n < 2: 

581 return 2 

582 if n < 7: 

583 return {2: 3, 3: 5, 4: 5, 5: 7, 6: 7}[n] 

584 if n <= sieve._list[-2]: 

585 l, u = sieve.search(n) 

586 if l == u: 

587 return sieve[u + 1] 

588 else: 

589 return sieve[u] 

590 nn = 6*(n//6) 

591 if nn == n: 

592 n += 1 

593 if isprime(n): 

594 return n 

595 n += 4 

596 elif n - nn == 5: 

597 n += 2 

598 if isprime(n): 

599 return n 

600 n += 4 

601 else: 

602 n = nn + 5 

603 while 1: 

604 if isprime(n): 

605 return n 

606 n += 2 

607 if isprime(n): 

608 return n 

609 n += 4 

610 

611 

612def prevprime(n): 

613 """ Return the largest prime smaller than n. 

614 

615 Notes 

616 ===== 

617 

618 Potential primes are located at 6*j +/- 1. This 

619 property is used during searching. 

620 

621 >>> from sympy import prevprime 

622 >>> [(i, prevprime(i)) for i in range(10, 15)] 

623 [(10, 7), (11, 7), (12, 11), (13, 11), (14, 13)] 

624 

625 See Also 

626 ======== 

627 

628 nextprime : Return the ith prime greater than n 

629 primerange : Generates all primes in a given range 

630 """ 

631 n = _as_int_ceiling(n) 

632 if n < 3: 

633 raise ValueError("no preceding primes") 

634 if n < 8: 

635 return {3: 2, 4: 3, 5: 3, 6: 5, 7: 5}[n] 

636 if n <= sieve._list[-1]: 

637 l, u = sieve.search(n) 

638 if l == u: 

639 return sieve[l-1] 

640 else: 

641 return sieve[l] 

642 nn = 6*(n//6) 

643 if n - nn <= 1: 

644 n = nn - 1 

645 if isprime(n): 

646 return n 

647 n -= 4 

648 else: 

649 n = nn + 1 

650 while 1: 

651 if isprime(n): 

652 return n 

653 n -= 2 

654 if isprime(n): 

655 return n 

656 n -= 4 

657 

658 

659def primerange(a, b=None): 

660 """ Generate a list of all prime numbers in the range [2, a), 

661 or [a, b). 

662 

663 If the range exists in the default sieve, the values will 

664 be returned from there; otherwise values will be returned 

665 but will not modify the sieve. 

666 

667 Examples 

668 ======== 

669 

670 >>> from sympy import primerange, prime 

671 

672 All primes less than 19: 

673 

674 >>> list(primerange(19)) 

675 [2, 3, 5, 7, 11, 13, 17] 

676 

677 All primes greater than or equal to 7 and less than 19: 

678 

679 >>> list(primerange(7, 19)) 

680 [7, 11, 13, 17] 

681 

682 All primes through the 10th prime 

683 

684 >>> list(primerange(prime(10) + 1)) 

685 [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] 

686 

687 The Sieve method, primerange, is generally faster but it will 

688 occupy more memory as the sieve stores values. The default 

689 instance of Sieve, named sieve, can be used: 

690 

691 >>> from sympy import sieve 

692 >>> list(sieve.primerange(1, 30)) 

693 [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] 

694 

695 Notes 

696 ===== 

697 

698 Some famous conjectures about the occurrence of primes in a given 

699 range are [1]: 

700 

701 - Twin primes: though often not, the following will give 2 primes 

702 an infinite number of times: 

703 primerange(6*n - 1, 6*n + 2) 

704 - Legendre's: the following always yields at least one prime 

705 primerange(n**2, (n+1)**2+1) 

706 - Bertrand's (proven): there is always a prime in the range 

707 primerange(n, 2*n) 

708 - Brocard's: there are at least four primes in the range 

709 primerange(prime(n)**2, prime(n+1)**2) 

710 

711 The average gap between primes is log(n) [2]; the gap between 

712 primes can be arbitrarily large since sequences of composite 

713 numbers are arbitrarily large, e.g. the numbers in the sequence 

714 n! + 2, n! + 3 ... n! + n are all composite. 

715 

716 See Also 

717 ======== 

718 

719 prime : Return the nth prime 

720 nextprime : Return the ith prime greater than n 

721 prevprime : Return the largest prime smaller than n 

722 randprime : Returns a random prime in a given range 

723 primorial : Returns the product of primes based on condition 

724 Sieve.primerange : return range from already computed primes 

725 or extend the sieve to contain the requested 

726 range. 

727 

728 References 

729 ========== 

730 

731 .. [1] https://en.wikipedia.org/wiki/Prime_number 

732 .. [2] https://primes.utm.edu/notes/gaps.html 

733 """ 

734 if b is None: 

735 a, b = 2, a 

736 if a >= b: 

737 return 

738 # if we already have the range, return it 

739 if b <= sieve._list[-1]: 

740 yield from sieve.primerange(a, b) 

741 return 

742 # otherwise compute, without storing, the desired range. 

743 

744 a = _as_int_ceiling(a) - 1 

745 b = _as_int_ceiling(b) 

746 while 1: 

747 a = nextprime(a) 

748 if a < b: 

749 yield a 

750 else: 

751 return 

752 

753 

754def randprime(a, b): 

755 """ Return a random prime number in the range [a, b). 

756 

757 Bertrand's postulate assures that 

758 randprime(a, 2*a) will always succeed for a > 1. 

759 

760 Examples 

761 ======== 

762 

763 >>> from sympy import randprime, isprime 

764 >>> randprime(1, 30) #doctest: +SKIP 

765 13 

766 >>> isprime(randprime(1, 30)) 

767 True 

768 

769 See Also 

770 ======== 

771 

772 primerange : Generate all primes in a given range 

773 

774 References 

775 ========== 

776 

777 .. [1] https://en.wikipedia.org/wiki/Bertrand's_postulate 

778 

779 """ 

780 if a >= b: 

781 return 

782 a, b = map(int, (a, b)) 

783 n = random.randint(a - 1, b) 

784 p = nextprime(n) 

785 if p >= b: 

786 p = prevprime(b) 

787 if p < a: 

788 raise ValueError("no primes exist in the specified range") 

789 return p 

790 

791 

792def primorial(n, nth=True): 

793 """ 

794 Returns the product of the first n primes (default) or 

795 the primes less than or equal to n (when ``nth=False``). 

796 

797 Examples 

798 ======== 

799 

800 >>> from sympy.ntheory.generate import primorial, primerange 

801 >>> from sympy import factorint, Mul, primefactors, sqrt 

802 >>> primorial(4) # the first 4 primes are 2, 3, 5, 7 

803 210 

804 >>> primorial(4, nth=False) # primes <= 4 are 2 and 3 

805 6 

806 >>> primorial(1) 

807 2 

808 >>> primorial(1, nth=False) 

809 1 

810 >>> primorial(sqrt(101), nth=False) 

811 210 

812 

813 One can argue that the primes are infinite since if you take 

814 a set of primes and multiply them together (e.g. the primorial) and 

815 then add or subtract 1, the result cannot be divided by any of the 

816 original factors, hence either 1 or more new primes must divide this 

817 product of primes. 

818 

819 In this case, the number itself is a new prime: 

820 

821 >>> factorint(primorial(4) + 1) 

822 {211: 1} 

823 

824 In this case two new primes are the factors: 

825 

826 >>> factorint(primorial(4) - 1) 

827 {11: 1, 19: 1} 

828 

829 Here, some primes smaller and larger than the primes multiplied together 

830 are obtained: 

831 

832 >>> p = list(primerange(10, 20)) 

833 >>> sorted(set(primefactors(Mul(*p) + 1)).difference(set(p))) 

834 [2, 5, 31, 149] 

835 

836 See Also 

837 ======== 

838 

839 primerange : Generate all primes in a given range 

840 

841 """ 

842 if nth: 

843 n = as_int(n) 

844 else: 

845 n = int(n) 

846 if n < 1: 

847 raise ValueError("primorial argument must be >= 1") 

848 p = 1 

849 if nth: 

850 for i in range(1, n + 1): 

851 p *= prime(i) 

852 else: 

853 for i in primerange(2, n + 1): 

854 p *= i 

855 return p 

856 

857 

858def cycle_length(f, x0, nmax=None, values=False): 

859 """For a given iterated sequence, return a generator that gives 

860 the length of the iterated cycle (lambda) and the length of terms 

861 before the cycle begins (mu); if ``values`` is True then the 

862 terms of the sequence will be returned instead. The sequence is 

863 started with value ``x0``. 

864 

865 Note: more than the first lambda + mu terms may be returned and this 

866 is the cost of cycle detection with Brent's method; there are, however, 

867 generally less terms calculated than would have been calculated if the 

868 proper ending point were determined, e.g. by using Floyd's method. 

869 

870 >>> from sympy.ntheory.generate import cycle_length 

871 

872 This will yield successive values of i <-- func(i): 

873 

874 >>> def iter(func, i): 

875 ... while 1: 

876 ... ii = func(i) 

877 ... yield ii 

878 ... i = ii 

879 ... 

880 

881 A function is defined: 

882 

883 >>> func = lambda i: (i**2 + 1) % 51 

884 

885 and given a seed of 4 and the mu and lambda terms calculated: 

886 

887 >>> next(cycle_length(func, 4)) 

888 (6, 2) 

889 

890 We can see what is meant by looking at the output: 

891 

892 >>> n = cycle_length(func, 4, values=True) 

893 >>> list(ni for ni in n) 

894 [17, 35, 2, 5, 26, 14, 44, 50, 2, 5, 26, 14] 

895 

896 There are 6 repeating values after the first 2. 

897 

898 If a sequence is suspected of being longer than you might wish, ``nmax`` 

899 can be used to exit early (and mu will be returned as None): 

900 

901 >>> next(cycle_length(func, 4, nmax = 4)) 

902 (4, None) 

903 >>> [ni for ni in cycle_length(func, 4, nmax = 4, values=True)] 

904 [17, 35, 2, 5] 

905 

906 Code modified from: 

907 https://en.wikipedia.org/wiki/Cycle_detection. 

908 """ 

909 

910 nmax = int(nmax or 0) 

911 

912 # main phase: search successive powers of two 

913 power = lam = 1 

914 tortoise, hare = x0, f(x0) # f(x0) is the element/node next to x0. 

915 i = 0 

916 while tortoise != hare and (not nmax or i < nmax): 

917 i += 1 

918 if power == lam: # time to start a new power of two? 

919 tortoise = hare 

920 power *= 2 

921 lam = 0 

922 if values: 

923 yield hare 

924 hare = f(hare) 

925 lam += 1 

926 if nmax and i == nmax: 

927 if values: 

928 return 

929 else: 

930 yield nmax, None 

931 return 

932 if not values: 

933 # Find the position of the first repetition of length lambda 

934 mu = 0 

935 tortoise = hare = x0 

936 for i in range(lam): 

937 hare = f(hare) 

938 while tortoise != hare: 

939 tortoise = f(tortoise) 

940 hare = f(hare) 

941 mu += 1 

942 if mu: 

943 mu -= 1 

944 yield lam, mu 

945 

946 

947def composite(nth): 

948 """ Return the nth composite number, with the composite numbers indexed as 

949 composite(1) = 4, composite(2) = 6, etc.... 

950 

951 Examples 

952 ======== 

953 

954 >>> from sympy import composite 

955 >>> composite(36) 

956 52 

957 >>> composite(1) 

958 4 

959 >>> composite(17737) 

960 20000 

961 

962 See Also 

963 ======== 

964 

965 sympy.ntheory.primetest.isprime : Test if n is prime 

966 primerange : Generate all primes in a given range 

967 primepi : Return the number of primes less than or equal to n 

968 prime : Return the nth prime 

969 compositepi : Return the number of positive composite numbers less than or equal to n 

970 """ 

971 n = as_int(nth) 

972 if n < 1: 

973 raise ValueError("nth must be a positive integer; composite(1) == 4") 

974 composite_arr = [4, 6, 8, 9, 10, 12, 14, 15, 16, 18] 

975 if n <= 10: 

976 return composite_arr[n - 1] 

977 

978 a, b = 4, sieve._list[-1] 

979 if n <= b - primepi(b) - 1: 

980 while a < b - 1: 

981 mid = (a + b) >> 1 

982 if mid - primepi(mid) - 1 > n: