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.

About these ads

Pages: 1 2

5 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

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 627 other followers

%d bloggers like this: