"""
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
