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)

About these ads

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 -- )
        2dup .header
        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),
    

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

%d bloggers like this: