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.

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 609 other followers

%d bloggers like this: