[Python] How to efficiently compute the exact value of prime counting function for n?

I want to know the exact count of prime numbers no greater than n, without finding all these prime numbers first.

I know of a very efficient method to find all primes no greater than n, using Sieve of Eratosthenes with Wheel Factorization, I use 2310-based wheel to skip all multiples of 2, 3, 5, 7, 11, so that I only need to process 480 numbers out of every 2310 numbers, skipping 61/77 of all candidates, and I only mark composites that are 2310-coprime and I only sieve up to square root, and I mark all multiples of 2, 3, 5, 7, 11 exactly once, so that my implementation is much much faster than naive Sieve of Eratosthenes found in online tutorials.

But even with my fast Sieve of Eratosthenes with Wheel Factorization computing all primes up to n for large n is slow. I have instead implemented one way to get the exact count of primes no greater than n without generating all the primes first.

I have implemented Meissel-Lehmer algorithm, and I based my implementation on this answer, I didn't implement the method from the other answer because it was way too complex.

I didn't use the recursive combinatorial phi(n, k) function because that is way too inefficient, instead, phi(n, k) means count of all numbers no greater than n coprime to all the first k primes, coprime means two numbers have no common divisors greater than 1, in other words phi(n, k) means the count of numbers that are not multiples of the first k prime numbers, and it can be efficiently calculated using my optimized prime sieve.

First, the combinatorial phi function implementation:

from bisect import bisect

class Wheel_Sieve:
def __init__(self):
self.primes = [2, 3, 5, 7, 11]
self.composites = {}
self.spoke = 0
self.coprime = 13
self.count = 5

def prime_range(self, limit: int) -> list:
self.sieve_by_value(limit)
return self.primes[: bisect(self.primes, limit)]

def prime_count(self, limit: int) -> int:
self.sieve_by_value(limit)
return bisect(self.primes, limit)

def first_primes(self, count: int) -> list:
self.sieve_by_count(count)
return self.primes[:count]

def sieve_by_count(self, limit: int) -> None:
while self.count < limit:
self.do_sieving(self.coprime)
self.coprime += WHEEL5[self.spoke]
self.spoke += 1
self.spoke -= 480 * (self.spoke == 480)

def sieve_by_value(self, limit: int) -> None:
while self.coprime <= limit:
self.do_sieving(self.coprime)
self.coprime += WHEEL5[self.spoke]
self.spoke += 1
self.spoke -= 480 * (self.spoke == 480)

def do_sieving(self, coprime: int) -> None:
duo = self.composites.pop(coprime, None)
if duo is None:
self.primes.append(coprime)
self.composites[coprime * coprime] = (self.spoke, coprime)
self.count += 1
else:
self.mark_composite(duo)

def mark_composite(self, duo):
j, prime = duo
multiple = self.coprime + prime * WHEEL5[j]
j += 1
j -= 480 * (j == 480)
while multiple in self.composites:
multiple += prime * WHEEL5[j]
j += 1
j -= 480 * (j == 480)
self.composites[multiple] = (j, prime)


SIEVE = Wheel_Sieve()

PHI_MEMO = {}


def phi_memo(n: int, k: int) -> int:
SIEVE.first_primes(k)
if result := PHI_MEMO.get((n, k)):
return result

if not k:
return n

if not n:
return 0

j = k - 1
result = phi_memo(n, j) - phi_memo(n // SIEVE.primes[j], j)
PHI_MEMO[(n, k)] = result
return result


And my efficient algorithm:

import numba as nb
import numpy as np


# fmt: off
FOREWHEEL5 = np.array((
(4 , 2 ), (9 , 6 ), (25 , 30 ), (35 , 30 ), (49 , 210 ), (77 , 210 ),
(91 , 210 ), (119 , 210 ), (133 , 210 ), (161 , 210 ), (203 , 210 ), (217 , 210 ),
(121 , 2310), (143 , 2310), (187 , 2310), (209 , 2310), (253 , 2310), (319 , 2310),
(341 , 2310), (407 , 2310), (451 , 2310), (473 , 2310), (517 , 2310), (583 , 2310),
(649 , 2310), (671 , 2310), (737 , 2310), (781 , 2310), (803 , 2310), (869 , 2310),
(913 , 2310), (979 , 2310), (1067, 2310), (1111, 2310), (1133, 2310), (1177, 2310),
(1199, 2310), (1243, 2310), (1331, 2310), (1397, 2310), (1441, 2310), (1507, 2310),
(1529, 2310), (1573, 2310), (1639, 2310), (1661, 2310), (1727, 2310), (1793, 2310),
(1837, 2310), (1859, 2310), (1903, 2310), (1969, 2310), (1991, 2310), (2057, 2310),
(2101, 2310), (2123, 2310), (2167, 2310), (2189, 2310), (2299, 2310), (2321, 2310)
), dtype=np.uint64)
WHEEL5 = np.array([
4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 6 , 6 , 2 , 6 , 4 , 2 , 6 , 4 , 6 , 8 , 4 , 2 , 4 , 2 , 4 ,
14, 4 , 6 , 2 , 10, 2 , 6 , 6 , 4 , 2 , 4 , 6 , 2 , 10, 2 , 4 , 2 , 12, 10, 2 , 4 , 2 , 4 , 6 ,
2 , 6 , 4 , 6 , 6 , 6 , 2 , 6 , 4 , 2 , 6 , 4 , 6 , 8 , 4 , 2 , 4 , 6 , 8 , 6 , 10, 2 , 4 , 6 ,
2 , 6 , 6 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 6 , 10, 2 , 10, 2 , 4 , 2 , 4 , 6 , 8 , 4 , 2 , 4 ,
12, 2 , 6 , 4 , 2 , 6 , 4 , 6 , 12, 2 , 4 , 2 , 4 , 8 , 6 , 4 , 6 , 2 , 4 , 6 , 2 , 6 , 10, 2 ,
4 , 6 , 2 , 6 , 4 , 2 , 4 , 2 , 10, 2 , 10, 2 , 4 , 6 , 6 , 2 , 6 , 6 , 4 , 6 , 6 , 2 , 6 , 4 ,
2 , 6 , 4 , 6 , 8 , 4 , 2 , 6 , 4 , 8 , 6 , 4 , 6 , 2 , 4 , 6 , 8 , 6 , 4 , 2 , 10, 2 , 6 , 4 ,
2 , 4 , 2 , 10, 2 , 10, 2 , 4 , 2 , 4 , 8 , 6 , 4 , 2 , 4 , 6 , 6 , 2 , 6 , 4 , 8 , 4 , 6 , 8 ,
4 , 2 , 4 , 2 , 4 , 8 , 6 , 4 , 6 , 6 , 6 , 2 , 6 , 6 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 2 ,
10, 2 , 10, 2 , 6 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 6 , 6 , 8 , 4 , 2 , 6 , 10, 8 , 4 , 2 , 4 , 2 ,
4 , 8 , 10, 6 , 2 , 4 , 8 , 6 , 6 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 6 , 2 , 10, 2 , 10, 2 , 4 , 2 ,
4 , 6 , 2 , 6 , 4 , 2 , 4 , 6 , 6 , 2 , 6 , 6 , 6 , 4 , 6 , 8 , 4 , 2 , 4 , 2 , 4 , 8 , 6 , 4 ,
8 , 4 , 6 , 2 , 6 , 6 , 4 , 2 , 4 , 6 , 8 , 4 , 2 , 4 , 2 , 10, 2 , 10, 2 , 4 , 2 , 4 , 6 , 2 ,
10, 2 , 4 , 6 , 8 , 6 , 4 , 2 , 6 , 4 , 6 , 8 , 4 , 6 , 2 , 4 , 8 , 6 , 4 , 6 , 2 , 4 , 6 , 2 ,
6 , 6 , 4 , 6 , 6 , 2 , 6 , 6 , 4 , 2 , 10, 2 , 10, 2 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 10, 6 ,
2 , 6 , 4 , 2 , 6 , 4 , 6 , 8 , 4 , 2 , 4 , 2 , 12, 6 , 4 , 6 , 2 , 4 , 6 , 2 , 12, 4 , 2 , 4 ,
8 , 6 , 4 , 2 , 4 , 2 , 10, 2 , 10, 6 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 6 , 6 , 2 , 6 , 4 , 2 ,
10, 6 , 8 , 6 , 4 , 2 , 4 , 8 , 6 , 4 , 6 , 2 , 4 , 6 , 2 , 6 , 6 , 6 , 4 , 6 , 2 , 6 , 4 , 2 ,
4 , 2 , 10, 12, 2 , 4 , 2 , 10, 2 , 6 , 4 , 2 , 4 , 6 , 6 , 2 , 10, 2 , 6 , 4 , 14, 4 , 2 , 4 ,
2 , 4 , 8 , 6 , 4 , 6 , 2 , 4 , 6 , 2 , 6 , 6 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 12, 2 , 12
], dtype=np.uint64)
# fmt: on


@nb.njit(cache=True)
def phi(n: int, limit: int) -> int:
sieve = np.full(n + 1, True, dtype=np.uint8)
sieve[:2] = False
start, step = FOREWHEEL5[0]
i = 0
j = 1
while start <= n and i < 60:
sieve[start::step] = False
start, step = FOREWHEEL5[j]
j += j != 59
i += 1

count = 5
last = 11
k = 13
i = 0
while count < limit and k <= n:
if sieve[k]:
count += 1
j = i
multiple = k * k
while multiple <= n:
sieve[multiple] = False
multiple += k * WHEEL5[j]
j += 1
j -= 480 * (j == 480)

last = k

k += WHEEL5
i += 1
i -= 480 * (i == 480)

return np.sum(sieve[last + 1 :]) + 1


My code is a thousand times faster:

In [364]: phi(65536, 1024)
Out[364]: 5519

In [365]: phi_memo(65536, 1024)
Out[365]: np.uint64(5519)

In [366]: %timeit PHI_MEMO.clear(); phi_memo(65536, 1024)
88.3 ms ± 371 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [367]: %timeit phi(65536, 1024)
91.5 μs ± 384 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)


I have implemented prime counting function using my phi function, it is indeed more efficient than sieving all the primes, but not by very much:

@nb.njit(cache=True)
def prime_sieve(limit: int) -> np.ndarray:
is_prime = np.full(limit + 1, True, dtype=np.uint8)
is_prime[:2] = False
start, step = FOREWHEEL5[0]
i = 0
j = 1
while start <= limit and i < 60:
is_prime[start::step] = False
start, step = FOREWHEEL5[j]
j += j != 59
i += 1

top = limit**0.5
low = int(top)
top = low + (top > low)
k = 13
i = 0
while k <= top:
if is_prime[k]:
j = i
multiple = k * k
while multiple <= limit:
is_prime[multiple] = False
multiple += k * WHEEL5[j]
j += 1
j -= 480 * (j == 480)

k += WHEEL5
i += 1
i -= 480 * (i == 480)

return is_prime


@nb.njit(cache=True)
def prime_sieve_with_pi(limit: int) -> tuple[np.ndarray, np.ndarray]:
is_prime = prime_sieve(limit)
return np.flatnonzero(is_prime), np.cumsum(is_prime.astype(np.uint64))


@nb.njit(cache=True)
def get_primes(limit: int) -> np.ndarray:
return np.flatnonzero(prime_sieve(limit))


@nb.njit(cache=True)
def count_primes(limit: int) -> int:
return np.sum(prime_sieve(limit).astype(np.uint64))


def prime_count(limit: int) -> int:
sqrt = int(limit**0.5)
cbrt = int(limit ** (1 / 3))
primes, counts = prime_sieve_with_pi(limit // cbrt + 3)
low = counts[cbrt + 1]
high = counts[sqrt + 1]
return (
phi(limit, low)
+ low
- sum(counts[limit // primes] - i for i in range(low, high))
- 1
)


In [369]: prime_sieve_with_pi(10**8)
Out[369]:
(array([ 2, 3, 5, ..., 99999959, 99999971, 99999989]),
array([ 0, 0, 1, ..., 5761455, 5761455, 5761455],
dtype=uint64))

In [370]: prime_count(10**8)
Out[370]: np.uint64(5761455)

In [371]: %timeit prime_count(10**8)
316 ms ± 2.64 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [372]: %timeit prime_sieve_with_pi(10**8)
963 ms ± 8.19 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [373]: %timeit prime_sieve(10**8)
377 ms ± 2.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [374]: %timeit count_primes(10**8)
637 ms ± 5.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [375]: %timeit get_primes(10**8)
502 ms ± 2.77 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [376]: count_primes(10**8)
Out[376]: 5761455


How can I compute the exact value of prime counting function much faster than sieving?

Continue reading...