## Factor Tables

### May 8, 2012

Sieving for least prime factors is a simple variant of the normal Sieve of Eratosthenes:

(define (lpf-sieve n)
(let* ((len (quotient (- n 1) 2))
(lpf (make-vector len 1)))
(let loop ((i 0) (p 3))
(cond ((< n (* p p)) lpf)
((= (vector-ref lpf i) 1)
(do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
((<= len j) (loop (+ i 1) (+ p 2)))
(when (= (vector-ref lpf j) 1)
(vector-set! lpf j p))))
(else (loop (+ i 1) (+ p 2)))))))

We sieve to a million:

(define lpf (lpf-sieve 1000000))

The factor table makes it easy to factor a number, using the same method as you would manually:

(define (factors n)
(let loop ((n n) (fs (list)))
(if (= n 1) (reverse fs)
(let ((f (if (even? n) 2 (vector-ref lpf (/ (- n 3) 2)))))
(if (= f 1) (reverse (cons n fs)) (loop (/ n f) (cons f fs)))))))

Displaying a table of least prime factors is just a matter of getting the loops right:

(define (show-table k n00)
(define (p->i p) (/ (- p 3) 2))
(display "  ")
(do ((j 0 (+ j 1))) ((= j k) (newline))
(show5 (+ j n00)))
(do ((i 0 (+ i 10))) ((= i 100))
(do ((ns '(1 3 7 9) (cdr ns))) ((null? ns))
(show2 (+ i (car ns)))
(do ((j 0 (+ j 1))) ((= j k))
(let ((p (+ (* 100 (+ j n00)) (+ i (car ns)))))
(if (= p 1) (display "     ")
(let ((f (vector-ref lpf (p->i p))))
(if (= f 1) (display "   --") (show5 f))))))
(newline))))

Scheme is noted for its lack of formatting capability. The functions that format the table are ugly, but simple:

(define (show2 num)
(if (< num 10) (display " "))
(display num))

(define (show5 num)
(if (< num 10) (display "    ")
(if (< num 100) (display "   ")
(if (< num 1000) (display "  ")
(if (< num 10000) (display " ")))))
(display num))

Here is the last page of a table of least prime factors to a million:

> (show-table 10 9990)
9990 9991 9992 9993 9994 9995 9996 9997 9998 9999
1   19   --    3  181   13    3   29  107    3    7
3    3    7   89    3   83  191    3   37    7    3
7   --   17    3   --   37    3    7   43    3   --
9    3  127   17    3   31    7    3  223   --    3
11   13    3   41   17    3   23   --    3  487   11
13  347   79    3    7   17    3   --   11    3   19
17  859    3   --   11    3  281   17    3    7   --
19    7   11    3  641   19    3  157    7    3  991
21    3  191   --    3   53   --    3   --   17    3
23   --    3   31   13    3    7   --    3   11   17
27    3  229  541    3   11   53    3   --  139    3
29   --    3    7   --    3   --   37    3  167    7
31   11    7    3   --   --    3   --  599    3   --
33    3   --   --    3   --   19    3    7   23    3
37   13   29    3  233   --    3   59   47    3  113
39    3  277   --    3    7   41    3   13  283    3
41   71    3   61    7    3   --  149    3   67  577
43   --   23    3   19   73    3   47  653    3    7
47    7    3  383   31    3  193   11    3  467   13
49   --   --    3   13   11    3    7   --    3   29
51    3   73   11    3   --    7    3   71   37    3
53   11    3   29  163    3   --   --    3   --   --
57    3  821    7    3   19   13    3   11  109    3
59  107    3   37   --    3   11   29    3    7   --
61    7   31    3   11  503    3   13    7    3   --
63    3   11  223    3  601   --    3   --   --    3
67   --   13    3  911    7    3   --  173    3   31
69    3   --   --    3  151   83    3   --   13    3
71   83    3    7   --    3   19   --    3  127    7
73   17    7    3   23  257    3  479   --    3   13
77   19    3   17   --    3  199    7    3  691   11
79   29  199    3   17   13    3  163   11    3   --
81    3   --   23    3    7   11    3   31   73    3
83   --    3   59    7    3   13   --    3   --   --
87    3    7   --    3  107   79    3   17    7    3
89    7    3  331   --    3   73  137    3   11   19
91   --   19    3   97   --    3    7   13    3   17
93    3   13   41    3   11    7    3   43   71    3
97   11   97    3    7  499    3  953  433    3  757
99    3   --    7    3   --   --    3   19  307    3

You can run the program at http://programmingpraxis.codepad.org/a9eDAmhr; the calculations are right, but the formatting gets mangled.

There’s more you can do with sieving that just compute primes or least prime factors. For instance, we saw a sieve for twin primes in a recent exercise. And it is easy to use a sieve to compute a range of totients: initialize the elements of a sieve to the numbers that they represent, and for each multiple of each prime, divide the number in the sieve by the sieving prime.

(define (totients n)
(let ((tots (make-vector (+ n 1))))
(do ((i 0 (+ i 1))) ((< n i))
(vector-set! tots i i))
(do ((i 2 (+ i 1))) ((< n i) tots)
(when (= i (vector-ref tots i))
(vector-set! tots i (- i 1))
(do ((j (+ i i) (+ i j))) ((< n j))
(vector-set! tots j
(* (vector-ref tots j) (- 1 (/ i)))))))))

> (totients 100)
#101(0 1 1 2 2 4 2 6 4 6 4 10 4 12 6 8 8 16 6 18 8 12 10 22
8 20 12 18 12 28 8 30 16 20 16 24 12 36 18 24 16 40 12 42 20
24 22 46 16 42 20 32 24 52 18 40 24 36 28 58 16 60 30 36 32
48 20 66 32 44 24 70 24 72 36 40 36 60 24 78 32 54 40 82 24
64 42 56 40 88 24 72 44 60 46 72 32 96 42 60 40)

Pages: 1 2

### 5 Responses to “Factor Tables”

1. David said

Forth implementation:

1000000 constant LPF_SIZE

: sq ( x -- x^2 )  s" dup *" evaluate ; immediate

: create-table  ( compile: sz -- ; run: n -- addr )
create
here over cells allot
swap 0 DO 1 !+ LOOP drop
does>
swap 1- cells + ;

LPF_SIZE create-table lpf

: sieve-lpf ( -- )
2
BEGIN
dup lpf @ 1 = IF
dup sq
BEGIN dup LPF_SIZE <= WHILE
dup lpf @ 1 = IF
2dup lpf !
THEN
over +
REPEAT drop
THEN
1+
dup sq LPF_SIZE > UNTIL
drop ;

: .factors  ( n -- )
BEGIN  dup lpf @  dup 1 <>  WHILE
dup .
/
REPEAT drop
. ;

: .header ( n00 k -- )
cr ."   "
0 DO  dup 5 .r  1+  LOOP  drop cr ;

: show?  ( n -- ? )
dup 1 and 0<>  swap 5 mod 0<>  and ;

: .row  ( n00 k row -- )
dup 2 .r -rot
bounds DO
i 100 * over +  lpf @
dup 1 = IF
."    --" drop
ELSE
5 .r
THEN
LOOP drop cr ;

: .table  ( n00 k -- )
100 0 DO
i show? IF 2dup i .row THEN
LOOP 2drop ;

sieve-lpf

Testing:

9990 10 .table
9990 9991 9992 9993 9994 9995 9996 9997 9998 9999
1   19   --    3  181   13    3   29  107    3    7
3    3    7   89    3   83  191    3   37    7    3
7   --   17    3   --   37    3    7   43    3   --
9    3  127   17    3   31    7    3  223   --    3
11   13    3   41   17    3   23   --    3  487   11
13  347   79    3    7   17    3   --   11    3   19
17  859    3   --   11    3  281   17    3    7   --
19    7   11    3  641   19    3  157    7    3  991
21    3  191   --    3   53   --    3   --   17    3
23   --    3   31   13    3    7   --    3   11   17
27    3  229  541    3   11   53    3   --  139    3
29   --    3    7   --    3   --   37    3  167    7
31   11    7    3   --   --    3   --  599    3   --
33    3   --   --    3   --   19    3    7   23    3
37   13   29    3  233   --    3   59   47    3  113
39    3  277   --    3    7   41    3   13  283    3
41   71    3   61    7    3   --  149    3   67  577
43   --   23    3   19   73    3   47  653    3    7
47    7    3  383   31    3  193   11    3  467   13
49   --   --    3   13   11    3    7   --    3   29
51    3   73   11    3   --    7    3   71   37    3
53   11    3   29  163    3   --   --    3   --   --
57    3  821    7    3   19   13    3   11  109    3
59  107    3   37   --    3   11   29    3    7   --
61    7   31    3   11  503    3   13    7    3   --
63    3   11  223    3  601   --    3   --   --    3
67   --   13    3  911    7    3   --  173    3   31
69    3   --   --    3  151   83    3   --   13    3
71   83    3    7   --    3   19   --    3  127    7
73   17    7    3   23  257    3  479   --    3   13
77   19    3   17   --    3  199    7    3  691   11
79   29  199    3   17   13    3  163   11    3   --
81    3   --   23    3    7   11    3   31   73    3
83   --    3   59    7    3   13   --    3   --   --
87    3    7   --    3  107   79    3   17    7    3
89    7    3  331   --    3   73  137    3   11   19
91   --   19    3   97   --    3    7   13    3   17
93    3   13   41    3   11    7    3   43   71    3
97   11   97    3    7  499    3  953  433    3  757
99    3   --    7    3   --   --    3   19  307    3
ok
32767 .factors 7 31 151  ok

2. Python implementation:

Testing:

3. ror said

c#: http://pastebin.com/41KxNmpP
I didn’t notice that numbers divisible by 2 and 5 are omitted, and I replaced — with blanks

4. Paul said

Another Python implementation.

pagenr = 1000  # 1000 numbers on a page; 1: 0-1000
N = 1000
last = pagenr * N
lowprime = ["--"] * (last + 1)
for i in range(3, last + 1, 2):
if lowprime[i] == "--":
for j in xrange(min(i * i, last), last + 1, 2 * i):
if lowprime[j] == "--":
lowprime[j] = i
fs = lambda x: "{0:>6s}".format(str(x))
print fs(""),
for k in xrange((last - N) // 100, (last) // 100):
print fs(k),
for row in xrange(1, 100, 2):
if row % 5:
print
print fs(row),
for k in lowprime[row + last - N::100]:
print fs(k),