Source code for sympy.ntheory.generate

"""
Generating and counting primes.

"""
from __future__ import print_function, division

import random
from bisect import bisect
# Using arrays for sieving instead of lists greatly reduces
# memory consumption
from array import array as _array

from sympy import Function, S
from sympy.core.compatibility import as_int, range
from .primetest import isprime


def _azeros(n):
    return _array('l', [0]*n)


def _aset(*v):
    return _array('l', v)


def _arange(a, b):
    return _array('l', range(a, b))


[docs]class Sieve: """An infinite list of prime numbers, implemented as a dynamically growing sieve of Eratosthenes. When a lookup is requested involving an odd number that has not been sieved, the sieve is automatically extended up to that number. Examples ======== >>> from sympy import sieve >>> sieve._reset() # this line for doctest only >>> 25 in sieve False >>> sieve._list array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23]) """ # data shared (and updated) by all Sieve instances def __init__(self): self._n = 6 self._list = _aset(2, 3, 5, 7, 11, 13) # primes self._tlist = _aset(0, 1, 1, 2, 2, 4) # totient self._mlist = _aset(0, 1, -1, -1, 0, -1) # mobius assert all(len(i) == self._n for i in (self._list, self._tlist, self._mlist)) def __repr__(self): return ("<%s sieve (%i): %i, %i, %i, ... %i, %i\n" "%s sieve (%i): %i, %i, %i, ... %i, %i\n" "%s sieve (%i): %i, %i, %i, ... %i, %i>") % ( 'prime', len(self._list), self._list[0], self._list[1], self._list[2], self._list[-2], self._list[-1], 'totient', len(self._tlist), self._tlist[0], self._tlist[1], self._tlist[2], self._tlist[-2], self._tlist[-1], 'mobius', len(self._mlist), self._mlist[0], self._mlist[1], self._mlist[2], self._mlist[-2], self._mlist[-1]) def _reset(self, prime=None, totient=None, mobius=None): """Reset all caches (default). To reset one or more set the desired keyword to True.""" if all(i is None for i in (prime, totient, mobius)): prime = totient = mobius = True if prime: self._list = self._list[:self._n] if totient: self._tlist = self._tlist[:self._n] if mobius: self._mlist = self._mlist[:self._n]
[docs] def extend(self, n): """Grow the sieve to cover all primes <= n (a real number). Examples ======== >>> from sympy import sieve >>> sieve._reset() # this line for doctest only >>> sieve.extend(30) >>> sieve[10] == 29 True """ n = int(n) if n <= self._list[-1]: return # We need to sieve against all bases up to sqrt(n). # This is a recursive call that will do nothing if there are enough # known bases already. maxbase = int(n**0.5) + 1 self.extend(maxbase) # Create a new sieve starting from sqrt(n) begin = self._list[-1] + 1 newsieve = _arange(begin, n + 1) # Now eliminate all multiples of primes in [2, sqrt(n)] for p in self.primerange(2, maxbase): # Start counting at a multiple of p, offsetting # the index to account for the new sieve's base index startindex = (-begin) % p for i in range(startindex, len(newsieve), p): newsieve[i] = 0 # Merge the sieves self._list += _array('l', [x for x in newsieve if x])
[docs] def extend_to_no(self, i): """Extend to include the ith prime number. Parameters ========== i : integer Examples ======== >>> from sympy import sieve >>> sieve._reset() # this line for doctest only >>> sieve.extend_to_no(9) >>> sieve._list array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23]) Notes ===== The list is extended by 50% if it is too short, so it is likely that it will be longer than requested. """ i = as_int(i) while len(self._list) < i: self.extend(int(self._list[-1] * 1.5))
[docs] def primerange(self, a, b): """Generate all prime numbers in the range [a, b). Examples ======== >>> from sympy import sieve >>> print([i for i in sieve.primerange(7, 18)]) [7, 11, 13, 17] """ from sympy.functions.elementary.integers import ceiling # wrapping ceiling in as_int will raise an error if there was a problem # determining whether the expression was exactly an integer or not a = max(2, as_int(ceiling(a))) b = as_int(ceiling(b)) if a >= b: return self.extend(b) i = self.search(a)[1] maxi = len(self._list) + 1 while i < maxi: p = self._list[i - 1] if p < b: yield p i += 1 else: return
[docs] def totientrange(self, a, b): """Generate all totient numbers for the range [a, b). Examples ======== >>> from sympy import sieve >>> print([i for i in sieve.totientrange(7, 18)]) [6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16] """ from sympy.functions.elementary.integers import ceiling # wrapping ceiling in as_int will raise an error if there was a problem # determining whether the expression was exactly an integer or not a = max(1, as_int(ceiling(a))) b = as_int(ceiling(b)) n = len(self._tlist) if a >= b: return elif b <= n: for i in range(a, b): yield self._tlist[i] else: self._tlist += _arange(n, b) for i in range(1, n): ti = self._tlist[i] startindex = (n + i - 1) // i * i for j in range(startindex, b, i): self._tlist[j] -= ti if i >= a: yield ti for i in range(n, b): ti = self._tlist[i] for j in range(2 * i, b, i): self._tlist[j] -= ti if i >= a: yield ti
[docs] def mobiusrange(self, a, b): """Generate all mobius numbers for the range [a, b). Parameters ========== a : integer First number in range b : integer First number outside of range Examples ======== >>> from sympy import sieve >>> print([i for i in sieve.mobiusrange(7, 18)]) [-1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1] """ from sympy.functions.elementary.integers import ceiling # wrapping ceiling in as_int will raise an error if there was a problem # determining whether the expression was exactly an integer or not a = max(1, as_int(ceiling(a))) b = as_int(ceiling(b)) n = len(self._mlist) if a >= b: return elif b <= n: for i in range(a, b): yield self._mlist[i] else: self._mlist += _azeros(b - n) for i in range(1, n): mi = self._mlist[i] startindex = (n + i - 1) // i * i for j in range(startindex, b, i): self._mlist[j] -= mi if i >= a: yield mi for i in range(n, b): mi = self._mlist[i] for j in range(2 * i, b, i): self._mlist[j] -= mi if i >= a: yield mi
[docs] def search(self, n): """Return the indices i, j of the primes that bound n. If n is prime then i == j. Although n can be an expression, if ceiling cannot convert it to an integer then an n error will be raised. Examples ======== >>> from sympy import sieve >>> sieve.search(25) (9, 10) >>> sieve.search(23) (9, 9) """ from sympy.functions.elementary.integers import ceiling # wrapping ceiling in as_int will raise an error if there was a problem # determining whether the expression was exactly an integer or not test = as_int(ceiling(n)) n = as_int(n) if n < 2: raise ValueError("n should be >= 2 but got: %s" % n) if n > self._list[-1]: self.extend(n) b = bisect(self._list, n) if self._list[b - 1] == test: return b, b else: return b, b + 1
def __contains__(self, n): try: n = as_int(n) assert n >= 2 except (ValueError, AssertionError): return False if n % 2 == 0: return n == 2 a, b = self.search(n) return a == b def __getitem__(self, n): """Return the nth prime number""" if isinstance(n, slice): self.extend_to_no(n.stop) return self._list[n.start - 1:n.stop - 1:n.step] else: n = as_int(n) self.extend_to_no(n) return self._list[n - 1]
# Generate a global object for repeated use in trial division etc sieve = Sieve()
[docs]def prime(nth): """ Return the nth prime, with the primes indexed as prime(1) = 2, prime(2) = 3, etc.... The nth prime is approximately n*log(n). Logarithmic integral of x is a pretty nice approximation for number of primes <= x, i.e. li(x) ~ pi(x) In fact, for the numbers we are concerned about( x<1e11 ), li(x) - pi(x) < 50000 Also, li(x) > pi(x) can be safely assumed for the numbers which can be evaluated by this function. Here, we find the least integer m such that li(m) > n using binary search. Now pi(m-1) < li(m-1) <= n, We find pi(m - 1) using primepi function. Starting from m, we have to find n - pi(m-1) more primes. For the inputs this implementation can handle, we will have to test primality for at max about 10**5 numbers, to get our answer. Examples ======== >>> from sympy import prime >>> prime(10) 29 >>> prime(1) 2 >>> prime(100000) 1299709 See Also ======== sympy.ntheory.primetest.isprime : Test if n is prime primerange : Generate all primes in a given range primepi : Return the number of primes less than or equal to n References ========== .. [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 .. [2] https://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number .. [3] https://en.wikipedia.org/wiki/Skewes%27_number """ n = as_int(nth) if n < 1: raise ValueError("nth must be a positive integer; prime(1) == 2") if n <= len(sieve._list): return sieve[n] from sympy.functions.special.error_functions import li from sympy.functions.elementary.exponential import log a = 2 # Lower bound for binary search b = int(n*(log(n) + log(log(n)))) # Upper bound for the search. while a < b: mid = (a + b) >> 1 if li(mid) > n: b = mid else: a = mid + 1 n_primes = primepi(a - 1) while n_primes < n: if isprime(a): n_primes += 1 a += 1 return a - 1
[docs]class primepi(Function): """ Represents the prime counting function pi(n) = the number of prime numbers less than or equal to n. Algorithm Description: In sieve method, we remove all multiples of prime p except p itself. Let phi(i,j) be the number of integers 2 <= k <= i which remain after sieving from primes less than or equal to j. Clearly, pi(n) = phi(n, sqrt(n)) If j is not a prime, phi(i,j) = phi(i, j - 1) if j is a prime, We remove all numbers(except j) whose smallest prime factor is j. Let x= j*a be such a number, where 2 <= a<= i / j Now, after sieving from primes <= j - 1, a must remain (because x, and hence a has no prime factor <= j - 1) Clearly, there are phi(i / j, j - 1) such a which remain on sieving from primes <= j - 1 Now, if a is a prime less than equal to j - 1, x= j*a has smallest prime factor = a, and has already been removed(by sieving from a). So, we don't need to remove it again. (Note: there will be pi(j - 1) such x) Thus, number of x, that will be removed are: phi(i / j, j - 1) - phi(j - 1, j - 1) (Note that pi(j - 1) = phi(j - 1, j - 1)) => phi(i,j) = phi(i, j - 1) - phi(i / j, j - 1) + phi(j - 1, j - 1) So,following recursion is used and implemented as dp: phi(a, b) = phi(a, b - 1), if b is not a prime phi(a, b) = phi(a, b-1)-phi(a / b, b-1) + phi(b-1, b-1), if b is prime Clearly a is always of the form floor(n / k), which can take at most 2*sqrt(n) values. Two arrays arr1,arr2 are maintained arr1[i] = phi(i, j), arr2[i] = phi(n // i, j) Finally the answer is arr2[1] Examples ======== >>> from sympy import primepi >>> primepi(25) 9 See Also ======== sympy.ntheory.primetest.isprime : Test if n is prime primerange : Generate all primes in a given range prime : Return the nth prime """ @classmethod def eval(cls, n): if n is S.Infinity: return S.Infinity if n is S.NegativeInfinity: return S.Zero try: n = int(n) except TypeError: if n.is_real == False or n is S.NaN: raise ValueError("n must be real") return if n < 2: return S.Zero if n <= sieve._list[-1]: return S(sieve.search(n)[0]) lim = int(n ** 0.5) lim -= 1 lim = max(lim, 0) while lim * lim <= n: lim += 1 lim -= 1 arr1 = [0] * (lim + 1) arr2 = [0] * (lim + 1) for i in range(1, lim + 1): arr1[i] = i - 1 arr2[i] = n // i - 1 for i in range(2, lim + 1): # Presently, arr1[k]=phi(k,i - 1), # arr2[k] = phi(n // k,i - 1) if arr1[i] == arr1[i - 1]: continue p = arr1[i - 1] for j in range(1, min(n // (i * i), lim) + 1): st = i * j if st <= lim: arr2[j] -= arr2[st] - p else: arr2[j] -= arr1[n // st] - p lim2 = min(lim, i * i - 1) for j in range(lim, lim2, -1): arr1[j] -= arr1[j // i] - p return S(arr2[1])
[docs]def nextprime(n, ith=1): """ Return the ith prime greater than n. i must be an integer. Notes ===== Potential primes are located at 6*j +/- 1. This property is used during searching. >>> from sympy import nextprime >>> [(i, nextprime(i)) for i in range(10, 15)] [(10, 11), (11, 13), (12, 13), (13, 17), (14, 17)] >>> nextprime(2, ith=2) # the 2nd prime after 2 5 See Also ======== prevprime : Return the largest prime smaller than n primerange : Generate all primes in a given range """ n = int(n) i = as_int(ith) if i > 1: pr = n j = 1 while 1: pr = nextprime(pr) j += 1 if j > i: break return pr if n < 2: return 2 if n < 7: return {2: 3, 3: 5, 4: 5, 5: 7, 6: 7}[n] if n <= sieve._list[-2]: l, u = sieve.search(n) if l == u: return sieve[u + 1] else: return sieve[u] nn = 6*(n//6) if nn == n: n += 1 if isprime(n): return n n += 4 elif n - nn == 5: n += 2 if isprime(n): return n n += 4 else: n = nn + 5 while 1: if isprime(n): return n n += 2 if isprime(n): return n n += 4
[docs]def prevprime(n): """ Return the largest prime smaller than n. Notes ===== Potential primes are located at 6*j +/- 1. This property is used during searching. >>> from sympy import prevprime >>> [(i, prevprime(i)) for i in range(10, 15)] [(10, 7), (11, 7), (12, 11), (13, 11), (14, 13)] See Also ======== nextprime : Return the ith prime greater than n primerange : Generates all primes in a given range """ from sympy.functions.elementary.integers import ceiling # wrapping ceiling in as_int will raise an error if there was a problem # determining whether the expression was exactly an integer or not n = as_int(ceiling(n)) if n < 3: raise ValueError("no preceding primes") if n < 8: return {3: 2, 4: 3, 5: 3, 6: 5, 7: 5}[n] if n <= sieve._list[-1]: l, u = sieve.search(n) if l == u: return sieve[l-1] else: return sieve[l] nn = 6*(n//6) if n - nn <= 1: n = nn - 1 if isprime(n): return n n -= 4 else: n = nn + 1 while 1: if isprime(n): return n n -= 2 if isprime(n): return n n -= 4
[docs]def primerange(a, b): """ Generate a list of all prime numbers in the range [a, b). If the range exists in the default sieve, the values will be returned from there; otherwise values will be returned but will not modify the sieve. Examples ======== >>> from sympy import primerange, sieve >>> print([i for i in primerange(1, 30)]) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] The Sieve method, primerange, is generally faster but it will occupy more memory as the sieve stores values. The default instance of Sieve, named sieve, can be used: >>> list(sieve.primerange(1, 30)) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29] Notes ===== Some famous conjectures about the occurrence of primes in a given range are [1]: - Twin primes: though often not, the following will give 2 primes an infinite number of times: primerange(6*n - 1, 6*n + 2) - Legendre's: the following always yields at least one prime primerange(n**2, (n+1)**2+1) - Bertrand's (proven): there is always a prime in the range primerange(n, 2*n) - Brocard's: there are at least four primes in the range primerange(prime(n)**2, prime(n+1)**2) The average gap between primes is log(n) [2]; the gap between primes can be arbitrarily large since sequences of composite numbers are arbitrarily large, e.g. the numbers in the sequence n! + 2, n! + 3 ... n! + n are all composite. See Also ======== nextprime : Return the ith prime greater than n prevprime : Return the largest prime smaller than n randprime : Returns a random prime in a given range primorial : Returns the product of primes based on condition Sieve.primerange : return range from already computed primes or extend the sieve to contain the requested range. References ========== .. [1] https://en.wikipedia.org/wiki/Prime_number .. [2] http://primes.utm.edu/notes/gaps.html """ from sympy.functions.elementary.integers import ceiling if a >= b: return # if we already have the range, return it if b <= sieve._list[-1]: for i in sieve.primerange(a, b): yield i return # otherwise compute, without storing, the desired range. # wrapping ceiling in as_int will raise an error if there was a problem # determining whether the expression was exactly an integer or not a = as_int(ceiling(a)) - 1 b = as_int(ceiling(b)) while 1: a = nextprime(a) if a < b: yield a else: return
[docs]def randprime(a, b): """ Return a random prime number in the range [a, b). Bertrand's postulate assures that randprime(a, 2*a) will always succeed for a > 1. Examples ======== >>> from sympy import randprime, isprime >>> randprime(1, 30) #doctest: +SKIP 13 >>> isprime(randprime(1, 30)) True See Also ======== primerange : Generate all primes in a given range References ========== .. [1] https://en.wikipedia.org/wiki/Bertrand's_postulate """ if a >= b: return a, b = map(int, (a, b)) n = random.randint(a - 1, b) p = nextprime(n) if p >= b: p = prevprime(b) if p < a: raise ValueError("no primes exist in the specified range") return p
[docs]def primorial(n, nth=True): """ Returns the product of the first n primes (default) or the primes less than or equal to n (when ``nth=False``). Examples ======== >>> from sympy.ntheory.generate import primorial, randprime, primerange >>> from sympy import factorint, Mul, primefactors, sqrt >>> primorial(4) # the first 4 primes are 2, 3, 5, 7 210 >>> primorial(4, nth=False) # primes <= 4 are 2 and 3 6 >>> primorial(1) 2 >>> primorial(1, nth=False) 1 >>> primorial(sqrt(101), nth=False) 210 One can argue that the primes are infinite since if you take a set of primes and multiply them together (e.g. the primorial) and then add or subtract 1, the result cannot be divided by any of the original factors, hence either 1 or more new primes must divide this product of primes. In this case, the number itself is a new prime: >>> factorint(primorial(4) + 1) {211: 1} In this case two new primes are the factors: >>> factorint(primorial(4) - 1) {11: 1, 19: 1} Here, some primes smaller and larger than the primes multiplied together are obtained: >>> p = list(primerange(10, 20)) >>> sorted(set(primefactors(Mul(*p) + 1)).difference(set(p))) [2, 5, 31, 149] See Also ======== primerange : Generate all primes in a given range """ if nth: n = as_int(n) else: n = int(n) if n < 1: raise ValueError("primorial argument must be >= 1") p = 1 if nth: for i in range(1, n + 1): p *= prime(i) else: for i in primerange(2, n + 1): p *= i return p
[docs]def cycle_length(f, x0, nmax=None, values=False): """For a given iterated sequence, return a generator that gives the length of the iterated cycle (lambda) and the length of terms before the cycle begins (mu); if ``values`` is True then the terms of the sequence will be returned instead. The sequence is started with value ``x0``. Note: more than the first lambda + mu terms may be returned and this is the cost of cycle detection with Brent's method; there are, however, generally less terms calculated than would have been calculated if the proper ending point were determined, e.g. by using Floyd's method. >>> from sympy.ntheory.generate import cycle_length This will yield successive values of i <-- func(i): >>> def iter(func, i): ... while 1: ... ii = func(i) ... yield ii ... i = ii ... A function is defined: >>> func = lambda i: (i**2 + 1) % 51 and given a seed of 4 and the mu and lambda terms calculated: >>> next(cycle_length(func, 4)) (6, 2) We can see what is meant by looking at the output: >>> n = cycle_length(func, 4, values=True) >>> list(ni for ni in n) [17, 35, 2, 5, 26, 14, 44, 50, 2, 5, 26, 14] There are 6 repeating values after the first 2. If a sequence is suspected of being longer than you might wish, ``nmax`` can be used to exit early (and mu will be returned as None): >>> next(cycle_length(func, 4, nmax = 4)) (4, None) >>> [ni for ni in cycle_length(func, 4, nmax = 4, values=True)] [17, 35, 2, 5] Code modified from: https://en.wikipedia.org/wiki/Cycle_detection. """ nmax = int(nmax or 0) # main phase: search successive powers of two power = lam = 1 tortoise, hare = x0, f(x0) # f(x0) is the element/node next to x0. i = 0 while tortoise != hare and (not nmax or i < nmax): i += 1 if power == lam: # time to start a new power of two? tortoise = hare power *= 2 lam = 0 if values: yield hare hare = f(hare) lam += 1 if nmax and i == nmax: if values: return else: yield nmax, None return if not values: # Find the position of the first repetition of length lambda mu = 0 tortoise = hare = x0 for i in range(lam): hare = f(hare) while tortoise != hare: tortoise = f(tortoise) hare = f(hare) mu += 1 if mu: mu -= 1 yield lam, mu
[docs]def composite(nth): """ Return the nth composite number, with the composite numbers indexed as composite(1) = 4, composite(2) = 6, etc.... Examples ======== >>> from sympy import composite >>> composite(36) 52 >>> composite(1) 4 >>> composite(17737) 20000 See Also ======== sympy.ntheory.primetest.isprime : Test if n is prime primerange : Generate all primes in a given range primepi : Return the number of primes less than or equal to n prime : Return the nth prime compositepi : Return the number of positive composite numbers less than or equal to n """ n = as_int(nth) if n < 1: raise ValueError("nth must be a positive integer; composite(1) == 4") composite_arr = [4, 6, 8, 9, 10, 12, 14, 15, 16, 18] if n <= 10: return composite_arr[n - 1] a, b = 4, sieve._list[-1] if n <= b - primepi(b) - 1: while a < b - 1: mid = (a + b) >> 1 if mid - primepi(mid) - 1 > n: b = mid else: a = mid if isprime(a): a -= 1 return a from sympy.functions.special.error_functions import li from sympy.functions.elementary.exponential import log a = 4 # Lower bound for binary search b = int(n*(log(n) + log(log(n)))) # Upper bound for the search. while a < b: mid = (a + b) >> 1 if mid - li(mid) - 1 > n: b = mid else: a = mid + 1 n_composites = a - primepi(a) - 1 while n_composites > n: if not isprime(a): n_composites -= 1 a -= 1 if isprime(a): a -= 1 return a
[docs]def compositepi(n): """ Return the number of positive composite numbers less than or equal to n. The first positive composite is 4, i.e. compositepi(4) = 1. Examples ======== >>> from sympy import compositepi >>> compositepi(25) 15 >>> compositepi(1000) 831 See Also ======== sympy.ntheory.primetest.isprime : Test if n is prime primerange : Generate all primes in a given range prime : Return the nth prime primepi : Return the number of primes less than or equal to n composite : Return the nth composite number """ n = int(n) if n < 4: return 0 return n - primepi(n) - 1