## 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.

Pages: 1 2

### 6 Responses to “Wheel Factorization”

1. […] Praxis – Wheel Factorization By Remco Niemeijer In today’s Programming Praxis problem we’re supposed to factor numbers both the naive way and using […]

2. Remco Niemeijer said

```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 600851475143
```
3. Remco Niemeijer said

Thanks 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 600851475143
```
4. There 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 [] 2
```
5. David said

My 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. ;
```
6. Ethan Lee said

``` 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] } ```

```} ```