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)
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 -- ) 2dup .header 100 0 DO i show? IF 2dup i .row THEN LOOP 2drop ; sieve-lpfTesting:
Python implementation:
Testing:
Python implementation: http://pastebin.com/CFP3S2HV
Testing: http://pastebin.com/T1Q0W8pP
c#: http://pastebin.com/41KxNmpP
I didn’t notice that numbers divisible by 2 and 5 are omitted, and I replaced — with blanks
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),