Wheel Factorisation

Here is a description of a Python implementation of Wheel Factorisation to find prime numbers.  The goal is to see how quickly π(n), the number of primes ≤ n, can be calculated for values up to π(1011) using  wheel factorisation, comparing various wheel sizes.

How wheel factorisation works

Wheel Factorisation arranges numbers to be sieved in the spokes of a wheel.  The number of spokes, W, is the product of the first few primes,2, 3, …, pw , so W = 2 x 3 x …x pw,  for example W = 2 x 3 x 5 = 30. Each spoke has an initial number m (1≤m≤W) and contains numbers of the form m + jW (j=0,1,2,…). In other words, each spoke contains the numbers congruent to m modulo W.

Spokes whose initial number has a common factor with W can be eliminated from the sieve, as numbers in those spokes (other than the primes 2,.., pw themselves) are all composites.

The spokes that need to be considered for sieving are therefore those whose initial number has no common factor with W. The following diagram shows these spokes for W = 2x3x5 = 30, the spokes being rearranged as rows in an array. There are 8 spokes, with initial numbers 1, 7, 11, 13, 17, 19, 23, 29. Each row of the array represents one congruence class modulo W. (The blue row at the top is not a spoke, but a counter for the columns)

This array is sieved in a similar way to the Sieve of Eratosthenes, by eliminating in turn multiples of each prime.

Multiples of a prime p > pw need to be eliminated from each row (congruence class), so we need to find a starting point in each congruence class. One way to do this is to multiply p by the initial number of each congruence class. The resulting numbers are guaranteed to be in different congruence classes, as the following analysis shows:

Suppose  a and b are in different congruence classes, i.e.  a ≠b (mod W)
⇒ (a-b) ≠0 (mod W)
⇒ (a-b)p ≠0 (mod W),  (as p and W are relatively prime)
⇒ ap ≠bp (mod W)

In our example of W=30,  multiplying p by 1, 7, 11, 13, 17, 19, 23, 29 gives numbers p, 7p, 11p, 13p, 17p, 19p, 23p, 29p.

But starting the sieve elimination with these numbers would eliminate p itself, so instead of multiplying p by 1, multiplying it by (1+W) instead will produce the next multiple of p in the same congruence class.

The diagram below shows 7p, 11p, 13p, 17p, 19p, 23p, 29p, 31p, for  p=41. Yellow numbers are numbers by which p=41 is multiplied. Red numbers are the resulting products.

Having found the initial multiple of p in each congruence class, we can then eliminate multiples of p in each equivalence class by eliminating each p-th entry in the equivalence class, starting with the entries highlighted in red.

But we can improve on this strategy. Multiples of p that are smaller than p2 need not be eliminated, as they will already have been eliminated when sieving smaller primes. So the starting point for each congruence class can be the first multiple of p that is ≥ p2. These starting points can be found by multiplying p by the first numbers in each congruence class that are ≥ p.

The following figure illustrates this technique for p=41. Yellow numbers are numbers by which p=41 is multiplied. Red numbers are the resulting products.

Python implementation

The Python implementation of this algorithm uses the bitarray extension to efficiently store a Boolean value indicating whether a number is prime. The algorithm returns a dictionary of bitarrays, with the values of the bitarrays defining whether or not a particular number is prime.

import bitarray
import fractions

def WheelOfPrimes(n, Baseprimes):

    W=1
    for p in Baseprimes:
      W*=p
    Hub = []

    if n<W*W+1:
        return [] # we're not going to bother for n n:
            spokes[i] = bitarray.bitarray(spokeLen-1)
        else:
            spokes[i] = bitarray.bitarray(spokeLen)
        spokes[i].setall(1)

    p=Hub[1]
    i,c = 1,0
    while p*p         if spokes[Hub[i]][c] :
          j,cc = i,c
          for m in range(H):
            xp_div, xp_mod =divmod((Hub[j]+cc*W)*p,W)
            spokes[xp_mod][xp_div::p]=0
            j=(j+1)%H
            if j==0:
                cc+=1

        # find the next potential prime in the wheel
        i=(i+1)%H
        if i==0:
           c+=1
        p = Hub[i]+c*W

    spokes[1][0]=0

    return spokes

The performance of the implementation is tested by calculating π(n), the number of primes ≤ n.  The values of π(n) for powers of 10 can be verified against published lists, for example the website How Many Primes Are There?

Example code to run the algorithm is:

from time import time
start_time = time()
base="[2,3,5,7]"
testnum=1000*1000*1000
w = WheelOfPrimes(testnum, base)
c=len(base)
for i in w:
  c+=w[i].count()
print testnum,',"',  base,'",', c,",", round(time() - start_time,3)

The following graph shows the time taken to calculate to calculate π(n) per billion for various bases and values of n between 500,000,000 and 500,000,000,000.

For larger values of n, the best bases are [2,3,5,7] and [2,3,5,7,11], so extending the limit of n to 1011 for these two bases, the absolute time to calculate π(n) looks like:

 

Here is a table of these values:

[2, 3, 5, 7, 11]

[2, 3, 5, 7]

5 x 108

6.0

3.0

1 x 109

8.4

4.9

1 x 1010

68.6

60.2

5 x 1010

344.5

338.3

1 x 1011

718.3

727.8

The best time to calculate π(1011) is 718.3s (11 min 58.3s), for a wheel size of 2 x 3 x 5 x 7 x 11 = 2310.

1011 is about the largest value for which I can use this algorithm on a 4GB laptop running Windows 7.

The algorithm lends itself to multithreading, which might substantially reduce the time taken, but alas not under Windows 7. If anyone reading this wants to experiment with multithreading the algorithm, I would be interested in the results.