2

I have already implemented Sieve of Eratosthenes and another version with Wheel factorization, both in NumPy. But I don't think I have done it right, because the version with wheel factorization is actually slower.

I have already made several optimizations to the sieve. First all even numbers except 2 are composites (multiples of 2), so I only check odd numbers, the number of candidates is halved. Second all prime numbers except 2 and 3 are of the form 6k+1 or 6k-1, because if the modulo 6 remainder is 3 then the number must be a multiple of 3.

So starting from 5, only numbers with remainder of 1 or 5 need to be checked, the increment between each iteration is 6, thus the number of candidates is only a sixth of the original.

Third the start of the multiple for each prime number found is its square, this ensures this multiple cannot be previously flipped, and the step of the multiple increment is 2*n, because even multiples are already eliminated, so only odd multiples remain.

This is the version without Wheel factorization:

import numpy as np

def prime_sieve(n: int) -> np.ndarray:
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = primes[4::2] = primes[9::6] = False
    for a in range(6, int(n**0.5) + 1, 6):
        for b in (-1, 1):
            i = a + b
            if primes[i]:
                primes[i * i :: 2 * i] = False

    return np.where(primes)[0]

It is highly inefficient, because it does way too many flipping operations, the same number can be flipped multiple times, and the flipping doesn't skip all the numbers that aren't of the forms 6k+1 and 6k-1 that are already eliminated, thus a lot of computational redundancy.

I have implemented another version with Wheel factorization, using information found on Wikipedia (I wrote the code completely by myself, there isn't pseudocode for Wheel sieve):

import numpy as np
from itertools import cycle

def prime_wheel_sieve(n: int) -> np.ndarray:
    wheel = cycle([4, 2, 4, 2, 4, 6, 2, 6])
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = primes[4::2] = False
    primes[9::6] = primes[15::10] = False
    k = 7
    while (square := k * k) <= n:
        if primes[k]:
            primes[square::2*k] = False
        k += next(wheel)
    return np.where(primes)[0]

But it suffers the same drawbacks as the one without Wheel, and is less performant:

In [412]: %timeit prime_sieve(65536)
226 µs ± 2.97 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [413]: %timeit prime_wheel_sieve(65536)
236 µs ± 24 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

How can I properly implement Sieve of Eratosthenes with Wheel factorization in NumPy?

Ξένη Γήινος
  • 2,181
  • 1
  • 9
  • 35
  • 2
    You're not really making use of Numpy as you should: just because you use Numpy arrays as your underlying data type doesn't make your program fast. You have to use Numpy functions/methods instead of Python loops and operations. Otherwise it'll be really slow. Whether it's actually feasible to implement what you want without Python loops is another story. – David May 29 '23 at 14:38
  • One random tip: you can massively improve your performance by creating `primes` with a `dtype` of `np.bool8`: `primes = np.ones(n + 1, dtype=np.bool8)`. – David May 29 '23 at 21:12
  • @David What is the advantage of using that `dtype` over normal `bool`? – jared Jun 12 '23 at 04:43
  • 1
    @jared I made the suggestion when no dtype was specified, so the default dtype (`np.float32`, I believe) was used. The question has since been updated. Relative to bool, there is no advantage. In fact, bool is better as I just noticed that `np.bool8` has been deprecated. So yes, use bool. – David Jun 12 '23 at 10:08
  • [Here](https://stackoverflow.com/questions/65969131/time-complexity-of-sieve-algorithm) you can find a wheel sieve version of size 30, to have a speed advantage n must be > 10^8. [Here](https://stackoverflow.com/questions/72756206/double-segmented-primes-sieve) instead you can find a segmented version of the same sieve of size 6. – user140242 Jun 22 '23 at 12:46

0 Answers0