Wheel Factorization
May 8, 2009
Here is a simple function to find factors by trial division:
(define (td-factors n)
(let loop ((n n) (x 2) (fs '()))
(cond ((< n (* x x)) (reverse (cons n fs)))
((zero? (modulo n x)) (loop (/ n x) x (cons x fs)))
(else (loop n (+ x 1) fs)))))
The list of totatives of a number is calculated by counting down from the number to zero and including in the output all those numbers that are coprime to the input:
(define (totatives n)
(let loop ((i n) (ts '()))
(cond ((zero? i) ts)
((= (gcd i n) 1)
(loop (- i 1) (cons i ts)))
(else (loop (- i 1) ts)))))
It is easier to work with the gaps between spokes than with the absolute positions of the spokes. For instance, the 2,3,5-wheel with spokes at 1, 7, 11, 13, 17, 19, 23 and 29 can be viewed as a wheel with gaps 6, 4, 2, 4, 2, 4, 6, and 2 (which wraps back around to 1). Diffs makes this calculation:
(define (diffs xs)
(let loop ((x (car xs)) (xs (cdr xs)) (ds '()))
(if (null? xs) (reverse ds)
(loop (car xs) (cdr xs) (cons (- (car xs) x) ds)))))
The wheel function builds a wheel, including its starting tail:
(define (wheel n)
(let* ((ps (primes n))
(t (apply * (cdr (reverse ps))))
(ts (totatives t))
(ds (diffs ts)))
(append (diffs ps)
(cycle (append (cdr ds)
(list 2)
(list (car ds)))))))
Wheel uses last-pair and cycle, and primes from a previous exercise:
(define (last-pair xs)
(if (null? (cdr xs)) xs
(last-pair (cdr xs))))
(define (cycle xs)
(set-cdr! (last-pair xs) xs) xs)
Finally, wheel-factors builds a 2,3,5,7-wheel and loops in the same manner as trial division. The call (wheel 11) builds a 2,3,5,7-wheel; 11 is the next prime after 7, and is needed to compute the starting point of the wheel:
(define wheel-factors
(let ((w (wheel 11)))
(lambda (n)
(let loop ((n n) (i 2) (fs '()) (w w))
(cond ((< n (* i i)) (reverse (cons n fs)))
((zero? (modulo n i))
(loop (quotient n i) i (cons i fs) w))
(else (loop n (+ i (car w)) fs (cdr w))))))))
Here we see the two factorization functions in action:
> (td-factors 600851475143)
(71 839 1471 6857)
> (wheel-factors 600851475143)
(71 839 1471 6857)
This program can be run at http://codepad.org/3kwNA12j.
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]
}
}