Primes in a Wide Range
Problem
Count the primes p with lo ≤ p ≤ hi. The bounds can be as large as 10^12, far too large to sieve from 1, but the window itself is narrow: hi minus lo is at most 10^6. Values below 2 inside the window are not prime and simply do not count.
Examples
Input: lo = 10, hi = 30
Output: 6
Why: 11, 13, 17, 19, 23 and 29Input: lo = 1000000000000, hi = 1000000001000
Output: 37
Why: only the 1,001 values of the window are sieved, not the trillion below themInput: lo = 14, hi = 16
Output: 0
Why: edge case, a window with no primes in itHints
0 / 3
Testing each value in the window by trial division costs up to a million divisions per value. Think about which divisors can ever matter for a number no bigger than hi.
A composite number up to hi always has a prime factor no larger than the square root of hi, which is at most 10^6 here. Those small primes can be found with an ordinary sieve.
Sieve the small primes up to the square root of hi. Then keep one flag per value of the window, and for each small prime p cross out its multiples inside the window, starting from the larger of p squared and the first multiple of p that is at least lo. Count the flags that survive.
Solution
Every composite value up to hi has a prime factor no larger than the square root of hi, so an ordinary sieve up to that root supplies every prime that can ever cross something out. The window then gets its own array of flags, indexed by value minus lo, and each small prime marks its multiples inside the window, starting at p squared or at the first multiple of p not below lo, whichever is later, so no prime crosses itself out. Slice assignment on a bytearray does each run of crossings in one step. Time is O(sqrt(hi) log log hi plus (hi minus lo) log log hi) and space is O(sqrt(hi) plus hi minus lo).
from math import isqrt
def primes_between(lo, hi):
lo = max(lo, 2) # 0 and 1 are not prime
if lo > hi:
return 0
root = isqrt(hi)
small, base = bytearray([1]) * (root + 1), []
for p in range(2, root + 1): # ordinary sieve up to sqrt(hi)
if small[p]:
base.append(p)
small[p * p::p] = bytearray(len(small[p * p::p]))
alive = bytearray([1]) * (hi - lo + 1) # one flag per value in the window
for p in base:
start = max(p * p, (lo + p - 1) // p * p)
alive[start - lo::p] = bytearray(len(alive[start - lo::p]))
return sum(alive)
print(primes_between(10, 30)) # -> 6
print(primes_between(10**12, 10**12 + 1000)) # -> 37
print(primes_between(14, 16)) # -> 0
print(primes_between(1, 10)) # -> 4Stuck on the idea rather than the code? Sieve of Eratosthenes covers it.