

Wheel factorization is a method to find prime numbers. So that’s what this post is about.

So what is this wheel factorization? First, say all you had was division to find primes. You start at 2 on the number line with an empty box for collected primes. Divide the current number by all the numbers in the box – if you have remainders from every division, the number you tested is prime. Since there’s no numbers in the box, 2 is prime. The next value on the line is 3. 3 / 2 gives you 1 with a remainder of 1, so three’s a prime. The next value to check is 4: 4 / 2 is 2 with no remainder and 4/3 is 1 with remainder 1. Since you didn’t get a remainder for every division, 4 is not a prime. Now you go to 5.

def trial_v1(end): "finds primes by trial division. returns a list" primes = [] # a collection box for found primes for i in range(2, end + 1): # loop over all integers 2, 3, 4, 5, ... end + 1 if all(map(lambda x: i % x, primes)): # modulo divide i by every found prime primes.append(i) # if every division gave a remainder, i is prime return primes # timing to find all primes under 50,000 #>>> t0 = time.clock(); v1 = trial_v1(50000); time.clock() - t0 #18.48303964047911

But stepping across the number line one at a time is slow. Instead of starting at 2, start at 3 and take steps of 2. Now you’re checking only the odd values to see if they’re prime – You’ve cut out half of all the integers you need to check. You can picture this as a wheel having one spoke and a circumference length of 2 rolling down the number line.

def trial_v2(end): # faster than trial_v1; skips all even integers "finds primes by trial division. returns a list" primes = [2] # starting prime for i in range(3, end + 1, 2): # loop over odd integers 3, 5, 7, 9, ... end if all(map(lambda x: i % x, primes)): primes.append(i) return primes # timing to find all primes under 50,000 #>>> t0 = time.clock(); v2 = trial_v2(50000); time.clock() - t0 #9.210297254653703 #>>> v1 == v2 # verify with trial_v1 results #True

If you notice that there’s a step of four between 7 and 11 (as 9 is not prime), you can make a larger wheel with primes 2 and 3 (with a circumference of 2*3 = 6), start rolling it at 5 to take alternate steps of 2 and 4. With this wheel you’ve now cut out even more integers you need to check.

Now not every odd number is prime, the wheel just cuts some values out. Values that you can be sure are prime are ones that are smaller than the square of the largest prime used to make the wheel. To make a wheel, start out with a few primes – say 2, 3 and 5. The circumference of this wheel will be the product of these primes: 2*3*5 = 30. The wheel will start rolling at the next prime – in this case, 7.

def prod(seq, factor=1): "sequence -> product" for i in seq: factor *= i return factor def spoke_steps(primes): """returns list of steps to each wheel gap that start from the last value in primes""" stpt = primes.pop() # where the wheel starts (pop removes last item in list) whcf = prod(primes) # wheel's circumference # gaps are every number that are divisible by initial primes (composites) spokes = [] # locate where the spokes are for i in xrange(stpt, stpt + whcf + 1, 2): if not all(map(lambda x: i % x, primes)): continue # gap else: spokes.append(i) # spoke # find the steps needed to jump to each spoke # (beginning from the start of the wheel) steps = [] # last step returns to start of wheel for i, j in enumerate(spokes): if not i: continue steps.append(j - spokes[i - 1]) return steps

Now, you need to add spokes to the wheel. Find all the values between 7 and 37 (starting point and circumference + starting point) that are not divisible by 2, 3 or 5. In this case it’s the values 7, 11, 13, 17, 19, 23, 29, 31 and 37. These are where your spokes will go on the wheel – the differences between them are 4, 2, 4, 2, 4, 6, 2 and 6 – the steps that are taken to skip over non primes.

This 2 – 5 wheel cuts out 73.33% of all composites (non primes). This percentage is the efficiency of the wheel, which is found by dividing the number of gaps on the wheel (where the spokes aren’t) by the circumference of the wheel. Since there are 8 spokes, there’s 30-8 = 22 gaps and the efficiency is 22/30.

def trial_v3(initial_primes, end): "finds primes by using a wheel factor and trial division. returns a list" steps = spoke_steps(initial_primes[:]) primes = initial_primes[:] i = primes[-1] # starting point si, sn = 0, len(steps) - 1 # step index, number of steps while 1: if all(map(lambda x: i % x, primes)): primes.append(i) i += steps[si] # roll wheel to next integer si += 1 # increase step index if si > sn: si = 0 # rollover step index if i > end: break return primes # timing to find all primes under 50,000 #>>> init_p = [2, 3, 5] # initial primes to build wheel #>>> t0 = time.clock(); v3 = trial_v3(init_p, 50000); time.clock() - t0 #6.050525852309605 #>>> v2 == v3 #True # more timings #>>> init_p = [2, 3, 5] # 2 - 3 wheel #>>> t0 = time.clock(); v3 = trial_v3(init_p, 50000); time.clock() - t0 #6.050525852309605 #>>> init_p = [2, 3, 5, 7] # 2 - 5 wheel #>>> t0 = time.clock(); v4 = trial_v3(init_p, 50000); time.clock() - t0 #4.885780981255678 # #>>> init_p = [2, 3, 5, 7, 11] # 2 - 7 wheel #>>> t0 = time.clock(); v5 = trial_v3(init_p, 50000); time.clock() - t0 #4.185585313178848 #>>> init_p = [2, 3, 5, 7, 11, 13] # 2 - 11 wheel #>>> t0 = time.clock(); v6 = trial_v3(init_p, 50000); time.clock() - t0 #3.848427562906977 #>>> init_p = [2, 3, 5, 7, 11, 13, 17] # 2 - 13 wheel #>>> t0 = time.clock(); v7 = trial_v3(init_p, 50000); time.clock() - t0 #3.579197986760107 #>>> t0 = time.clock(); v8 = trial_v3(init_p, 50000); time.clock() - t0 #3.7229256307713

You can make larger wheels, this one is a 2 – 7 wheel with circumference 210:

However, larger wheels quickly become impractical. A 2 – 23 wheel has a circumference of 9,699,690 and has 8,040,810 spokes. It’s best to pair a small wheel with another method of finding primes, like the sieve of Eratosthenes. Doing so will really get you movin’.