Refactoring a prime number generator generator

Discussion of Common Lisp

Refactoring a prime number generator generator

Postby Harleqin » Thu Jan 29, 2009 9:06 pm

I have written a little prime number generator generator (sic), but it is somehow ugly. I think that I need some help refactoring it, or someone pointing out what I could do more concise or efficient.

Description:

The generator returns a function which, whenever called, spits out the next prime number. For example, when it is first called, it returns 2, when it is called again, it returns 3, then 5 etc..

I implemented the sieve of Erastothenes by using a "wheel". This wheel is a list of lists or, to be more precise, it is an alist of numbers associated with their prime factors. Actually, the numbers are associated with a list of the double of each prime factor, since I only look at odd numbers (handling 2 as a special case at the beginning). The wheel always contains all numbers bigger than the current and for them, all prime factors (doubled) smaller than the current.

When the function is called, it first increments current, then checks it against the first number in the wheel. If it is equal, it is not prime, if not, then it is. When it is a prime, I begin to "strike out" all its multiples by putting it into the wheel. When it is not a prime, this part of the wheel is broken up by putting the prime factors on their respective next multiple; then the function recurs, until a prime is found.

Code: Select all
(defun make-prime-generator ()
  "returns a function with closure that gives all primes sequentially.
   Example: (let ((gimme-prime (make-prime-generator)))
              (funcall gimme-prime)  ; => 2
              (funcall gimme-prime)) ; => 3"
  (let ((current 1)
        ;; wheel is a list of lists. The sublists contain the
        ;; number and then all prime factors, sorted ascendingly.
        (wheel ()))
    (labels ((propagate-wheel ()
               (do ((wr (cdr wheel))
                    (wl wheel)          ; insert between wl and wr
                    (p (cdar wheel)))
                   ((null p)
                    (setf wheel (cdr wheel))
                    (next-prime))
                 (let ((strikeout (+ current (car p))))
                   (cond ((null wr)
                          (setf (cdr wl) (cons strikeout (car p))
                                wl (cdr wl)
                                p (cdr p)))
                         ((< strikeout (caar wr))
                          (setf (cdr wl) (cons (list strikeout (car p))
                                               wr)
                                wl (cdr wl)
                                p (cdr p)))
                         ((= strikeout (caar wr))
                          (setf (cdar wr) (merge 'list (cdar wr) (list (car p)) #'<))
                          (setf p (cdr p)))
                         ((> strikeout (caar wr))
                          (setf wr (cdr wr)
                                wl (cdr wl)))))))
             (put-prime-into-wheel (place number &optional (wl wheel))
               (cond ((null wheel)    ; outer wheel, when current is 3
                      (setf wheel (list (list place number))))
                     ((= (caar wl) place) ; inner
                      (setf (cdar wl) (merge 'list (cdar wl) (list number) #'<)))
                     ((null (cdr wl))
                      (setf (cdr wl) (list (list place number))))
                     ((< place (caadr wl))
                      (setf (cdr wl) (cons (list place number) (cdr wl))))
                     (t
                      (put-prime-into-wheel place number (cdr wl)))))
             (next-prime ()
               (case current
                 (1 (return-from next-prime (incf current)))
                 (2 (incf current))
                 (otherwise (incf current 2)))
               (if (and (car wheel)
                        (= (caar wheel) current))
                   ;; not a prime -> rotate wheel
                   (propagate-wheel)    ; will call next-prime again
                   ;; prime -> put it into the wheel, and return it
                   (progn
                     (put-prime-into-wheel (* 3 current) (* 2 current))
                     current))))
      #'next-prime)))


Is using an alist sensible here? Should I use a hash table instead?
"Just throw more hardware at it" is the root of all evil.
Svante
Harleqin
 
Posts: 71
Joined: Wed Dec 17, 2008 5:18 am
Location: Bonn, Germany

Re: Refactoring a prime number generator generator

Postby raffaelcavallaro » Mon Feb 02, 2009 3:10 pm

Code: Select all
? (defun make-prime-generator ()
   (let ((most-recent-prime 1))
     (lambda ()
       (labels ((prime-p (n)
                  (not (loop for i from 2 to (isqrt n)
                         when (integerp (/ n i))
                         do (return t)))))
         (loop for j from (1+ most-recent-prime)
           when (prime-p j)
           do (setf most-recent-prime j) (return j))))))
MAKE-PRIME-GENERATOR
? (let ((gimme-primes (make-prime-generator)))
    (dotimes (n 10) (print (funcall gimme-primes))))

2
3
5
7
11
13
17
19
23
29
NIL
? (let ((another-gimme-primes (make-prime-generator)))
    (dotimes (n 20) (print (funcall another-gimme-primes))))

2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
NIL
?
raffaelcavallaro
 
Posts: 5
Joined: Mon Feb 02, 2009 3:00 pm

Re: Refactoring a prime number generator generator

Postby raffaelcavallaro » Mon Feb 02, 2009 11:46 pm

this is a little clearer
Code: Select all
(defun make-prime-generator ()
    (let ((most-recent-prime 1))
      (lambda ()
        (flet ((prime-p (n)
                 (loop for i from 2 to (isqrt n)
                   when (integerp (/ n i))
                   return nil))))
        (loop for j from (1+ most-recent-prime)
          when (prime-p j)
          do (setf most-recent-prime j) (return j)))))
raffaelcavallaro
 
Posts: 5
Joined: Mon Feb 02, 2009 3:00 pm

Re: Refactoring a prime number generator generator

Postby Harleqin » Tue Feb 03, 2009 7:24 am

Sorry, but that doesn't work.

What does your prime-p return if it doesn't find an integer divisor? You put a closing parenthesis too many on line 7, too.

Even if it did work, the intended algorithm is much worse in runtime than the sieve of Erastothenes.

Alas, your first version works, and after testing both the primitive method and mine with the wheel through "time", I found that while the primitive method uses ten times more memory, they are almost even when it comes to processor cycles. This indicates that indeed, I have a bottleneck, and I strongly suspect that my alist is responsible. I have to traverse this list for each non-prime, after all. I will try the hash table instead.
"Just throw more hardware at it" is the root of all evil.
Svante
Harleqin
 
Posts: 71
Joined: Wed Dec 17, 2008 5:18 am
Location: Bonn, Germany

Re: Refactoring a prime number generator generator

Postby Harleqin » Tue Feb 03, 2009 8:53 am

OK, here is the hash version. It is much more concise.

Code: Select all
(defun make-prime-generator ()
  "returns a function with closure that gives all primes sequentially.
   Example: (let ((gimme-prime (make-prime-generator)))
              (funcall gimme-prime)  ; => 2
              (funcall gimme-prime)) ; => 3"
  (let ((current 1)
        ;; wheel is a hash-table of lists, number => prime-factors
        (wheel (make-hash-table)))
    (labels ((propagate-wheel ()
               (dolist (p (gethash current wheel))
                 (let ((strikeout (+ current p)))
                   (setf (gethash strikeout wheel)
                         (cons p (gethash strikeout wheel)))))
               (remhash current wheel)
               (next-prime))
             (put-prime-into-wheel (place number)
               (setf (gethash place wheel)
                     (cons number (gethash place wheel))))
             (next-prime ()
               (case current
                 (1 (return-from next-prime (incf current)))
                 (2 (incf current))
                 (otherwise (incf current 2)))
               (if (gethash current wheel)
                   ;; not a prime -> rotate wheel
                   (propagate-wheel)    ; will call next-prime again
                   ;; prime -> put it into the wheel, and return it
                   (progn
                     (put-prime-into-wheel (* 3 current) (* 2 current))
                     current))))
      #'next-prime)))


For comparison, the timings from summing the first 1000 primes:

Code: Select all
CL-USER> (time (test-pg (make-slow-prime-generator) 1000))
Evaluation took:
  0.044 seconds of real time
  0.036002 seconds of total run time (0.036002 user, 0.000000 system)
  81.82% CPU
  58,995,000 processor cycles
  1,265,664 bytes consed
3682913
CL-USER> (time (test-pg (make-alist-prime-generator) 1000))
Evaluation took:
  0.058 seconds of real time
  0.044003 seconds of total run time (0.044003 user, 0.000000 system)
  75.86% CPU
  76,906,674 processor cycles
  122,856 bytes consed
3682913
CL-USER> (time (test-pg (make-prime-generator) 1000))
Evaluation took:
  0.005 seconds of real time
  0.004000 seconds of total run time (0.004000 user, 0.000000 system)
  80.00% CPU
  6,249,442 processor cycles
  98,264 bytes consed
3682913


And the first 2000 primes:

Code: Select all
CL-USER> (time (test-pg (make-slow-prime-generator) 2000))
Evaluation took:
  0.114 seconds of real time
  0.100006 seconds of total run time (0.100006 user, 0.000000 system)
  87.72% CPU
  153,564,336 processor cycles
  3,735,544 bytes consed
16274627
CL-USER> (time (test-pg (make-alist-prime-generator) 2000))
Evaluation took:
  0.195 seconds of real time
  0.180012 seconds of total run time (0.180012 user, 0.000000 system)
  92.31% CPU
  261,008,527 processor cycles
  282,744 bytes consed
16274627
CL-USER> (time (test-pg (make-prime-generator) 2000))
Evaluation took:
  0.010 seconds of real time
  0.008000 seconds of total run time (0.008000 user, 0.000000 system)
  80.00% CPU
  14,036,611 processor cycles
  230,496 bytes consed
16274627


That will teach me to use the right data structure.
"Just throw more hardware at it" is the root of all evil.
Svante
Harleqin
 
Posts: 71
Joined: Wed Dec 17, 2008 5:18 am
Location: Bonn, Germany

Re: Refactoring a prime number generator generator

Postby raffaelcavallaro » Tue Feb 03, 2009 12:01 pm

Sorry, bad edit. This one works again:

Code: Select all
(defun make-prime-generator ()
    (let ((most-recent-prime 1))
      (lambda ()
        (flet ((prime-p (n)
                 (not (loop for i from 2 to (isqrt n)
                   when (integerp (/ n i))
                   return t))))
        (loop for j from (1+ most-recent-prime)
          when (prime-p j)
          do (setf most-recent-prime j) (return j))))))
raffaelcavallaro
 
Posts: 5
Joined: Mon Feb 02, 2009 3:00 pm


Return to Common Lisp

Who is online

Users browsing this forum: pjstirling and 4 guests