Ankeet Parikh

Prime sieve using NumPy

I came across a neat way to write sieves in Python using NumPy:

import numpy as np
def pi(N):
  b = np.zeros(N, dtype=bool)
  b[2:] = True
  for i in range(2, int(N**0.5 + 1)):
    if b[i]:
      b[i*i::i] = False
  return np.count_nonzero(b)
print(pi(10**2))

Note that all composite numbers less than $i^2$ have a prime factor less than $i$ so those will have been crossed out in previous iterations. The code uses the handy start:stop:step syntax, as well as scalar assignment to a slice. The runtime complexity of this algorithm is $O(n \log \log n )$. The constant factor appears to be pretty small, probably due to things like NumPy's C-implementation, an efficient memory representation, and various parallelisms. Running this on an Apple M1 Pro with 16GB memory, here are my profiling results: \[ \displaystyle \begin{array}{|c|c|c|} \hline N & \pi (N) & \text{time (sec)} \\ \hline 10^2 & 25 & 0.0000 \\ 10^3 & 168 & 0.0000 \\ 10^4 & 1229 & 0.0000 \\ 10^5 & 9592 & 0.0001 \\ 10^6 & 78498 & 0.0011 \\ 10^7 & 664579 & 0.0125 \\ 10^8 & 5761455 & 0.3574 \\ 10^9 & 50847534 & 5.5850 \\ 10^{10} & 455052511 & 143.7051 \\ \hline \end{array} \] There are much faster ways to count primes, but 6 seconds for $10^9$ isn't bad.