Python implementation of 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.

Using the optimum wheel size, π(1011) can be calculated in under 12 minutes.

Wheel Factorisation

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.




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 proves:

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.

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?

The following graph shows the time taken to calculate to calculate π(n)  for various bases and values of n.

	[2, 3, 5, 7, 11, 13]	[2, 3, 5, 7, 11]	[2, 3, 5, 7]	[2, 3, 5]	[2, 3]
1,000,000,000	52.8	8.0	5.2	6.2	9.4
5,000,000,000	152.9	34.9	29.6	34.5	42.4
10,000,000,000	182.4	64.7	61.4	76.4	91.3
50,000,000,000	723.4	329.8	321.4	419.0	518.3
100,000,000,000	1275.7	691.0	699.4	861.9

Time to calculate π(n) (seconds)

So the best time to calculate π(1011) is 691s (11 min 31s), for a wheel size of 2x3x5x7x11=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 multi threading, which might substantially reduce the time taken, but alas not under Windows 7. If anyone reading this wants to experiment with multi threading the algorithm, I would be interested in the results.

Python source for the algorithm is:

import bitarray

import fractions

def WheelOfPrimes(n, Baseprimes):

  

    W=1

    for p in Baseprimes:

      W*=p

    Hub = []

    if n

        return [] # we're not going to bother for n

        

    for i in range(1,W,2):

        if fractions.gcd(i,W)==1:

            Hub=Hub+[i]

    H=len(Hub)

    spokes = {}

    # spokes[i][j] represents then number  i + j*W

    # so the number m is represented by spokes[m%W][m//W]

    spokeLen = n//W+1

    for i in Hub:

        if i + (spokeLen-1)*W > 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 <=n :

        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

Example code to run WheelOfPrimes 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)