prime optimation
John Thingstad
Anyhow thought I would take anothe shot at it and came up with..
(defun slow-prime-p (number)
(cond
((member number '(0 1)) nil)
((= number 2) t)
((= (mod number 2) 0) nil)
(t
(loop for i from 3 to (sqrt number) by 2 do
(when (= (mod number i) 0)
(return-from slow-prime-p nil)))
t)))
(defun slow-count-primes< (number)
(loop for i to number counting (slow-prime-p i)))
(defun prime-p (number)
(declare (fixnum number)
(optimize (speed 3) (safety 0) #+lispworks (fixnum-safety 1)))
(cond
((<= number 10) (if (member number '(2 3 5 7)) t nil))
((zerop (rem number 2)) nil)
((zerop (rem number 3)) nil)
((zerop (rem number 5)) nil)
((zerop (rem number 7)) nil)
(t
(loop for i fixnum from 11 to (the fixnum (isqrt number)) by 2 do
(when (zerop (the fixnum (rem number i)))
(return-from prime-p nil)))
t)))
(defun count-primes< (number)
(declare (fixnum number))
(let ((count 0))
(declare (fixnum count))
(dotimes (i number)
(declare (fixnum i))
(when (prime-p i) (incf count)))
count))
CL-USER 19 > (time (slow-count-primes< 1000000))
Timing the evaluation of (SLOW-COUNT-PRIMES< 1000000)
User time = 10.453
System time = 0.000
Elapsed time = 10.500
Allocation = 12012260 bytes
0 Page faults
78498
CL-USER 20 > (time (count-primes< 1000000))
Timing the evaluation of (COUNT-PRIMES< 1000000)
User time = 0.953
System time = 0.000
Elapsed time = 0.953
Allocation = 524 bytes
0 Page faults
78498
Any better?
--
Using Opera's revolutionary e-mail client: http://www.opera.com/mail/
Harold Lee
[snip]
> CL-USER 20 > (time (count-primes< 1000000))
> Timing the evaluation of (COUNT-PRIMES< 1000000)
>
> User time = 0.953
> System time = 0.000
> Elapsed time = 0.953
> Allocation = 524 bytes
> 0 Page faults
> 78498
>
> Any better?
Here's what I use (SBCL 1.0, Linux 2.6, Pentium M 1.6 GHz). It uses (a
little) more memory in exchange for speed.
(defun primes (max &key (as-vector nil))
"Find the list of all primes p <= max."
(let ((prime-divisors (make-array (list max)
:element-type 'bit
:initial-element 1)))
;; All values default to prime (i.e. 1) in the prime-divisors
bitmap,
;; now we'll mark the non-primes 0 using a sieve of Eratosthenes.
(loop for i from 2 to (floor (sqrt max)) do
(loop for j from (* i 2) to max by i do
(setf (aref prime-divisors (1- j)) 0)))
(if as-vector (progn (setf (aref prime-divisors 0) 0)
prime-divisors)
;; Convert the results into a list of numbers.
(let ((nums (list 2)))
(loop for i from 3 to max by 2 do
(if (= 1 (aref prime-divisors (1- i))) (push i nums)))
(nreverse nums)))))
(defun count-primes< (n)
;; n >= 2
(1+ (loop with primes = (primes n :as-vector t)
for i from 3 to n by 2
sum (aref primes (1- i)))))
CL-USER> (time (count-primes< 1000000))
Evaluation took:
0.258 seconds of real time
0.254961 seconds of user run time
0.001 seconds of system run time
0 calls to %EVAL
0 page faults and
125,008 bytes consed.
78498
CL-USER>
Pillsy
[...]
> Here's what I use (SBCL 1.0, Linux 2.6, Pentium M 1.6 GHz). It uses (a
> little) more memory in exchange for speed.
> (defun primes (max &key (as-vector nil))
> "Find the list of all primes p <= max."
> (let ((prime-divisors (make-array (list max)
> :element-type 'bit
> :initial-element 1)))
> ;; All values default to prime (i.e. 1) in the prime-divisors
> :: bitmap, now we'll mark the non-primes 0 using a sieve of
> :: Eratosthenes.
> (loop for i from 2 to (floor (sqrt max)) do
You're doing more than twice as much work as you need to here. If i
reaches a number that has prime-divisor set to 0, you don't have to
eliminate all its multiples---they've already been eliminated. Using
this version speeds things up by more than a factor of 2:
(defun primes (max &key (as-vector nil))
"Find the list (or vector, if AS-VECTOR is true) of all
primes p <= MAX."
(let ((prime-divisors (make-array (list max)
:element-type 'bit
:initial-element 1)))
(loop
:for i :from 2 :to (isqrt max)
;; This test saves time!
:when (= 1 (aref prime-divisors (1- i)))
:do (loop
:for j :from (* i 2) :to max :by i
:do (setf (aref prime-divisors (1- j)) 0)))
(if as-vector
(progn (setf (aref prime-divisors 0) 0)
prime-divisors)
(cons 2
(loop
:for i :from 3 :to max :by 2
:when (= (aref prime-divisors (1- i)) 1)
:collect i)))))
If you call PRIMES repeatedly, it might be worthwhile to memoize it,
i.e., to keep the vector PRIME-DIVISORS around from one call to the
next.
Cheers,
Pillsy