I'm one of those who up to now hasn't been a fan of the extended loop
macro. Recently I decided that I might as well dig in and learn loop.
Then I thought, "I'll try one last time to wrap my mind around the
series package instead."
So I started reading the manual documents and started thinking about
trying to program the sieve algorithm using series. I've had a very
hard time of it. Has anyone done the sieve algorithm using series?
While mulling it over, I came up with the following `functional sieve
algorithm' which I thought was kind of nifty.
(defun make-factor-test (n)
"Return a function that identifies factors of n."
(declare (fixnum n))
(lambda (i)
(declare (fixnum i))
;; (format t "Checking if ~d is a factor of ~D~%" n i)
(zerop (mod i n))))
(defun functional-sieve (n)
(declare (fixnum n))
(do ((i 2 (1+ i)) ;; Start with 2.
(factor-test-list '())
(total-primes 0)
(last-divisor-needed (ceiling (sqrt n))))
((> i n) (format
t
"Total primes: ~D; length of factor test list: ~D~%"
total-primes (length factor-test-list)))
(declare (fixnum i total-primes last-divisor-needed))
(unless (some #'(lambda (factor-test) (funcall factor-test i))
factor-test-list)
;; Found a prime.
;; (format t "~D is prime.~%" i)
(when (< i last-divisor-needed)
(setf factor-test-list
;; Keep factor tests for smaller n in front.
(append factor-test-list (list (make-factor-test i)))))
(incf total-primes))))
I like this because it's simple. I think it's also very readable.
Unfortunately it's also slow (about a factor of 10) compared to the
sieve algorithm that uses a bit array and loop. It seems too simple
to optimize much but maybe there's something I missed.
One thing I tried was to make the `factor-test' functions work by
counting, so they wouldn't take an argument but would just count the
number of times they were called and return t every nth call. But
`some' is short-circuiting. :-)
--
Fred Gilham ······@csl.sri.com
The TMI accident was unique: it was the only multi-billion dollar
accident in history in which nobody was harmed.
-Howard C. Hayden, Professor Emeritus, U of Conn.
Fred Gilham <······@snapdragon.csl.sri.com> writes:
> Has anyone done the sieve algorithm using series?
I thought about this, and it's not obvious how to do it, as the sieve
uses the primes it's generating, so it's consuming its output. This
doesn't fit into the series model in any obvious way.
> While mulling it over, I came up with the following `functional sieve
> algorithm' which I thought was kind of nifty.
>
> [make a list of (zerop (mod i n)) tests and apply each]
Series are sort of functional, and my thoughts took a similar
direction at first. Since infinite series are lazy, I can just
repeatedly (SETF ints (CHOOSE-IF #'mod-test (SUBSERIES ints 1))), sort
of.
(defun func-sieve ()
(map-fn '(series integer)
#'(lambda (s) (collect-first s))
(scan-fn '(series integer)
#'(lambda () (scan-range :from 2))
#'(lambda (s)
(declare (type (series integer) s))
(choose-if #'(lambda (x)
(not (zerop (mod x (collect-first s)))))
(subseries s 1))))))
To get around the problem of consuming its own output, FUNC-SIEVE
generates a series of series (this is the SCAN-FN bit), representing
the intermediate stages of the computation: Each series starts from a
prime and has had the multiples of all lower primes removed.
CL-USER 1> (subseries (func-sieve) 0 100)
#Z(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89
97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179
181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271
277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379
383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479
487 491 499 503 509 521 523 541)
However, as you point out, this approach has the problem of using MOD
instead of +, so it's not a true sieve. It doesn't transform into
loops, either.
I think a true sieve has to have an array that it can stride over with
simple addition. That means using ALTER. Here's a straight-forward
translation of the usual array sieve. The strides are done by
(SCAN-RANGE :FROM (+ P P) :BY P).
(defun sieve (k)
(let ((primes (make-sequence `vector k)))
(alter (scan 'vector primes)
(scan-range))
(iterate ((p (choose (subseries (scan 'vector primes) 2 (isqrt k)))))
(alter (scan 'vector primes)
(mapping ((int (scan 'vector primes))
(m (mask (scan-range :from (+ p p) :by p))))
(if m nil int))))
(choose (subseries (scan 'vector primes) 2))))
This does get transformed into loops (apart from the last CHOOSE). It
could be optimized in the usual ways, omitting even numbers from the
array and using a bit array, but the idea is clearer this way.
--
Pekka P. Pirinen
If it's spam, it's a scam. Don't do business with net abusers.
> I think a true sieve has to have an array that it can stride over
> with simple addition. That means using ALTER. Here's a
> straight-forward translation of the usual array sieve. The strides
> are done by (SCAN-RANGE :FROM (+ P P) :BY P).
>
> (defun sieve (k)
> (let ((primes (make-sequence `vector k)))
> (alter (scan 'vector primes)
> (scan-range))
> (iterate ((p (choose (subseries (scan 'vector primes) 2 (isqrt k)))))
> (alter (scan 'vector primes)
> (mapping ((int (scan 'vector primes))
> (m (mask (scan-range :from (+ p p) :by p))))
> (if m nil int))))
> (choose (subseries (scan 'vector primes) 2))))
>
> This does get transformed into loops (apart from the last CHOOSE).
> It could be optimized in the usual ways, omitting even numbers from
> the array and using a bit array, but the idea is clearer this way.
Thanks for thinking about this. I'm impressed that you were able to
figure it out!
I think a barrier to using the series package is the lack of examples.
It never occurred to me to use alter, and your code makes it clearer
what it's for.
Using series is definitely brain-stretching, and that's one reason I
enjoy using Lisp.
--
Fred Gilham ······@csl.sri.com
The density of a textbook must be inversely proportional to the
density of the students using it. --- Dave Stringer-Calvert
>>>>> "Fred" == Fred Gilham <······@snapdragon.csl.sri.com> writes:
>> I think a true sieve has to have an array that it can stride over
>> with simple addition. That means using ALTER. Here's a
>> straight-forward translation of the usual array sieve. The strides
>> are done by (SCAN-RANGE :FROM (+ P P) :BY P).
>>
>> (defun sieve (k)
>> (let ((primes (make-sequence `vector k)))
>> (alter (scan 'vector primes)
>> (scan-range))
>> (iterate ((p (choose (subseries (scan 'vector primes) 2 (isqrt k)))))
>> (alter (scan 'vector primes)
>> (mapping ((int (scan 'vector primes))
>> (m (mask (scan-range :from (+ p p) :by p))))
>> (if m nil int))))
>> (choose (subseries (scan 'vector primes) 2))))
>>
>> This does get transformed into loops (apart from the last CHOOSE).
>> It could be optimized in the usual ways, omitting even numbers from
>> the array and using a bit array, but the idea is clearer this way.
Fred> Thanks for thinking about this. I'm impressed that you were able to
Fred> figure it out!
Fred> I think a barrier to using the series package is the lack of examples.
Fred> It never occurred to me to use alter, and your code makes it clearer
Fred> what it's for.
I'll try to add some examples to series. If other people have nice,
interesting examples, please send them to the series mailing list at
series.sourceforge.net. (You can also send them directly to me, but I
prefer the mailing list, lest I lose it by accident.)
Ray
Oops! ···············@globalgraphics.com (Pekka P. Pirinen) was seen spray-painting on a wall:
> Fred Gilham <······@snapdragon.csl.sri.com> writes:
>> Has anyone done the sieve algorithm using series?
>
> I thought about this, and it's not obvious how to do it, as the sieve
> uses the primes it's generating, so it's consuming its output. This
> doesn't fit into the series model in any obvious way.
Here's the implementation that I've seen previously:
(defun primes (max)
(series:choose-if
#'(lambda (n)
(series:collect-and
(series:mapping ((x (series:scan-range :from 3 :by 2
:upto (+ 1 (isqrt n)))))
(not (zerop (mod n x))))))
(series:scan-range :from 3 :by 2 :upto max)))
0]] (primes 500)
#Z(3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103
107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199
211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313
317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433
439 443 449 457 461 463 467 479 487 491 499)
There is, in effect, a "non-recursive" use of scan-range; the
implementation is not /ideal/ in that it is /not/ self-referencing.
It's still pretty quick :-) And I don't think the order of complexity
is fundamentally changed by the non-reflection.
--
(concatenate 'string "cbbrowne" ·@acm.org")
http://cbbrowne.com/info/lisp.html
"In computing, invariants are ephemeral." -- Alan Perlis