Wheel Factorization
May 8, 2009
Having discussed prime numbers in several previous exercises, we are now interested in the problem of factoring an integer n; for instance, the prime factors of 42 are 2, 3, and 7. A simple factoring method is to perform trial division by all the integers counting from 2 to the square root of n. Your first task is to write that function.
An easy optimization is to divide only by 2 and then by odd integers greater than 2, which saves half the work. A better optimization is to divide by 2, then 3, then 5, and thereafter to alternately add 2 and 4 to the trial divisors — 7, 11, 13, 17, 19, 23, and so on — since all prime numbers greater than 3 are of the form 6k±1 for some integer k.
It turns out that both those optimizations are special cases of a technique called wheel factorization. Consider a 2-wheel of circumference 2 rolling along a number line with a “spoke” at the number 1; if you start with the spoke at 3 on the number line, it will strike the number line at 5, then 7, and then every odd number after that. Or consider a 2,3-wheel of circumference 2×3=6 with spokes at the number 1 and 5; if you start with the 5-spoke at 5 on the number line, it will strike the number line at 7, 11, 13, 17, 19, 23, and so on. Or consider a 2,3,5-wheel of circumference 2×3×5=30 with spokes at 1, 7, 11, 13, 17, 19, 23 and 29 starting with the 29-spoke at 7. And so on: next is a 2,3,5,7-wheel, then a 2,3,5,7,11-wheel, and the sequence continues infinitely.
Wheel factorization works by performing trial division at each place where a spoke touches the number line. As the wheels grow larger, more and more of the trial divisors are prime, so less and less unnecessary work is done. Of course, there is a point of diminishing returns; when the wheel gets too large, it is just as much work to compute the wheel as to compute the list of primes, and costs just as much to store. But a small wheel is easy to compute, and not too big, and provides a simple optimization over naive trial division.
The spokes of the wheel are computed by looking for co-primes, which are those numbers for which the spoke has no factors in common with the circumference of the wheel; in other words, where the greatest common divisor of the spoke and the circumference is 1. For instance, a 2,3,5-wheel has a spoke at 17 because the greatest common divisor of 17 and 30 is 1, but no spoke at 18 because the greatest common divisor of 18 and 30 is 6. These numbers are called totatives; if you’re curious about the math behind them, ask your favorite search engine for information about Euler’s totient function.
It is easy to see this visually. Here is a list of the positive integers to 42, with primes highlighted:
1 2 3 4 5 6
7 8 9 10 11 12
13 14 15 16 17 18
19 20 21 22 23 24
25 26 27 28 29 30
31 32 33 34 35 36
37 38 39 40 41 42
After the first row, all the primes are in two columns, which correspond to the two spokes of a 2,3-wheel. If 853 were input to the 2,3-wheel factorization function, we would trial divide by 2, 3, 5, 7, 11, 13, 17, 19, 23, 25,and 29 before concluding that 853 was prime; note that 25 is not prime, but is relatively prime to the circumference of the wheel.
Your second task is to write a function that finds the factors of a given number using wheel factorization; you should compute and use a 2,3,5,7-wheel. What are the factors of 600851475143? When you are finished, you are welcome to read or run a suggested solution, or post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] Praxis – Wheel Factorization By Remco Niemeijer In today’s Programming Praxis problem we’re supposed to factor numbers both the naive way and using […]
My Haskell solution (see http://bonsaicode.wordpress.com/2009/05/08/programming-praxis-%E2%80%93-wheel-factorization/ for a version with comments):
import Data.List tdFactors :: Integer -> [Integer] -> [Integer] tdFactors n = nub . factor n . takeWhile (\x -> x * x < n) where factor _ [] = [] factor r (x:xs) | mod r x == 0 = x : factor (div r x) (x:xs) | otherwise = factor r xs factorNaive :: Integer -> [Integer] factorNaive n = tdFactors n [2..] spokes :: Integer -> [Bool] spokes c = map ((== 1) . gcd c) [1..c] wheel :: [Integer] -> [Integer] wheel ps = (ps ++) . drop 1 . map snd . filter fst $ zip (cycle . spokes $ product ps) [1..] factorWheel :: Integer -> [Integer] factorWheel n = tdFactors n $ wheel [2,3,5,7] main :: IO () main = do print $ factorNaive 600851475143 print $ factorWheel 600851475143Thanks to a suggestion from programmingpraxis I was able to speed up my program. Since tdFactors already more or less assumes a sorted list of integers as input we can stop looking as soon as x * x > r. The new code becomes:
tdFactors :: Integer -> [Integer] -> [Integer] tdFactors n = factor n where factor _ [] = [] factor r (x:xs) | x * x > r = [r] | mod r x == 0 = x : factor (div r x) (x:xs) | otherwise = factor r xs factorNaive :: Integer -> [Integer] factorNaive n = tdFactors n [2..] spokes :: Integer -> [Bool] spokes c = map ((== 1) . gcd c) [1..c] wheel :: [Integer] -> [Integer] wheel ps = (ps ++) . drop 1 . map snd . filter fst $ zip (cycle . spokes $ product ps) [1..] factorWheel :: Integer -> [Integer] factorWheel n = tdFactors n $ wheel [2,3,5,7] main :: IO () main = do print $ factorNaive 600851475143 print $ factorWheel 600851475143There isn’t much point in showing the code to compute the wheel; once it is calculated, it can be converted to a static (i.e., compile-time) value:
let wheel_factors = let wheel2357 = let rec w = 2 :: 4 :: 2 :: 4 :: 6 :: 2 :: 6 :: 4 :: 2 :: 4 :: 6 :: 6 :: 2 :: 6 :: 4 :: 2 :: 6 :: 4 :: 6 :: 8 :: 4 :: 2 :: 4 :: 2 :: 4 :: 8 :: 6 :: 4 :: 6 :: 2 :: 4 :: 6 :: 2 :: 6 :: 6 :: 4 :: 2 :: 4 :: 6 :: 2 :: 6 :: 4 :: 2 :: 4 :: 2 ::10 :: 2 ::10 :: w in 1 :: 2 :: 2 :: 4 :: w in let rec go (w :: ws as wheel) l p n = if p * p > n then List.rev (n :: l) else let d = n / p in if d * p == n then go wheel (p :: l) p d else go ws l (p + w) n in go wheel2357 [] 2My FORTH solution; modified brute force factorization we did in (I think a future) exercise to use the factor wheel. Cool thing is the defining word “cycle” to create a new data type — the cyclic constant array.
\ Quick word to create (non-cyclic, non-checked) constant array : const-array create does> swap cells + @ ; \ Create a cyclic data structure - constant array with \ automatic wrap-around. : cycle ( -- ) create here 0 , does> @+ >r swap r> mod cells + @ ; : end-cycle ( -- ) here over - cell / 1- swap ! ; cycle wheel2357 2 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 6 , 6 , 2 , 6 , 4 , 2 , 6 , 4 , 6 , 8 , 4 , 2 , 4 , 2 , 4 , 8 , 6 , 4 , 6 , 2 , 4 , 6 , 2 , 6 , 6 , 4 , 2 , 4 , 6 , 2 , 6 , 4 , 2 , 4 , 2 , 10 , 2 , 10 , end-cycle const-array start2357 1 , 2 , 2 , 4 , variable wheel_ndx : init_wheel ( -- ) -4 wheel_ndx ! ; : wheel_step ( -- n ) wheel_ndx @ dup 0< IF 4 + start2357 ELSE wheel2357 THEN 1 wheel_ndx +! ; : d/mod ( d n -- d%n d/n ) locals| n | \ FORTH does not have remainder function for double division \ so we multiply the quotient by the divisor and subtract... \ 2dup 1 n M*/ 2dup 2>r \ get quotient & save a copy n 1 M*/ d- drop \ compute remainder 2r> ; \ restore quotient : .factors ( d -- ) 2 locals| f | init_wheel BEGIN 2dup f dup m* d< not WHILE 2dup f d/mod rot 0= IF f . 2nip \ drop n; quotient on stack replaces it ELSE 2drop \ drop quotient f wheel_step + TO f THEN REPEAT d. ;public class PrimeCalc {
public static class PrimeWheel {
// http://primes.utm.edu/glossary/xpage/WheelFactorization.html
static final long[] PRIMES = { 2, 3, 5, 7 };
static final int WHEEL_FACTOR = 210; // 2 * 3 * 5 * 7;
static final long[] SIEVES = { 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149, 151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209, 211 };
public long calcFactor(long n) {
for (long p : PRIMES) {
if (n % p == 0) {
return p;
}
}
long lim = (long) StrictMath.sqrt(n);
for (long s = 0; s <= lim; s+= WHEEL_FACTOR) {
for (long sieve : SIEVES) {
long m = s + sieve;
if (n % m == 0) {
return m;
}
}
}
return n;
}
public boolean isPrime(long n) {
if (n < 2) {
return false;
}
return calcFactor(n) == n;
}
}
public static void main(String[] args) {
PrimeWheel pw = new PrimeWheel();
long num = 600851475143L;
java.util.List factors = new java.util.ArrayList();
while (true) {
long f = pw.calcFactor(num);
factors.add(f);
if (f == num) {
break;
} else {
num /= f;
}
}
System.out.println(factors);
// [71, 839, 1471, 6857]
}
}