Here's Óscar López's code, tweaked to build the list in the top-down manner:
(define (is-prime? n lop)
(define sqrtn (sqrt n))
(let loop ([lop lop])
(cond [(or (empty? lop) (< sqrtn (mcar lop))) true]
[(zero? (modulo n (mcar lop))) false]
[else (loop (mcdr lop))])))
(define (find-primes n)
(let* ([a (mcons 3 '())]
[b (mcons 2 a)])
(let loop ([p a] [i 5] [d 2] ; d = diff +2 +4 +2 ...
[c 2]) ; c = count of primes found
(cond [(> i n) c]
[(is-prime? i (mcdr a))
(set-mcdr! p (mcons i '()))
(loop (mcdr p) (+ i d) (- 6 d) (+ c 1))]
[else (loop p (+ i d) (- 6 d) c )]))))
Runs at about ~n1.25, empirically; compared to the original's ~n1.8..1.9, in the measured range, inside DrRacket (append is the culprit of that bad behaviour):
; (time (length (find-primes 100000))) ; with cons times in milliseconds
; 10K 156 ; 20K 437 ; 40K 1607 ; 80K 5241 ; 100K 7753 .... n^1.8-1.9-1.7 OP's
; 10K 62 ; 20K 109 ; 40K 421 ; 80K 1217 ; 100K 2293 .... n^1.8-1.9 Óscar's
; mcons:
(time (find-primes 2000000))
; 100K 47 ; 200K 172 ; 1M 1186 ; 2M 2839 ; 3M 4851 ; 4M 7036 .... n^1.25-1.32 this
; 9592 17984 78498 148933 216816 283146
It's still just a trial division though... :) The sieve of Eratosthenes will be much faster yet.
As for set-cdr!, it looks like we just have to use it, to emulate any lazy algorithm... otherwise, must have extendable arrays (lists of...), for the O(1) snoc/append1 operation; that's lots and lots of coding instead of just one primitive, frowned upon for some reason, unclear to me.
some more ramblings follow, concerning the history of the "sieve" in functional programming.
a historical footnote.
As naomik writes in the comments, Sussman/Abelson present a prime sieve in Structure and Interpretation of Computer Programs, Section 3.5.2. The code there looks similar to the previously known ones in FP circles, either that attributed to P. Quarendon in Henderson et. Morris's 1976 A Lazy Evaluator,
primeswrt[x;l] = if car[l] mod x=0 then primeswrt[x;cdr[l]] else
cons[car[l];primeswrt[x;cdr[l]]]
primes[l] = cons[car[l];primes[primeswrt[car[l];cdr[l]]]]
-- primes[integers[2]]
(somebody did implement the M-expressions, after all); or that of David Turner's,
from SASL 1976 manual. That I could only find, on the web, in its 1983 revision, with the rem code fully baked in.
But A.J.T.Davie 1992 (p. 136) gives it as
primes = sieve [2..]
sieve [p, ...nos...] = [p, ...sieve (remove (multsof p) nos)...]
(I'm writing this in a pattern-matching pseudocode, hopefully self-explanatory enough).
The question arises, what is multsof? Is it
multsof p = (n => rem n p == 0)
or perhaps
multsof p = [p, p+p..] -- ( or, even, [p*p, p*p+p..] )
? Both are possible; the choice fully determines the corresponding implementation of remove.
David Turner made the first choice, leading to remove being implemented as "a dual of filter". So have the authors of SICP.
But Abelson and Sussman, writing their code in the top-down manner as they often do
(though not for the sieve code), teach us a valuable lesson of delaying the implementation
(not just optimization) as much as possible.
The second choice would make the same code work as a bona fide sieve of Eratosthenes, instead of the trial division. Sieve of E. has much better complexity, so it should run noticeably faster for all but the smallest n.
Whatever the choice, the code itself has one more immediate concern, one pressing flaw -- it is too hasty with the removals, starts them way to early and ends up with a quadratic explosion in their numbers.
So the original, hasty
primes = sieve [2..]
sieve [p, ...nos... ] =
[p, ...sieve (remove (multsof p) nos)...]
should have been instead the slowly-yet-farther-progressing, postponing each removal until the right moment,
primes = [2, ...sieve [3..] primes... ]
sieve [...ns..., p*p, ...nos...] [p, ...ps...] = -- with non-linear
[...ns..., ...sieve (remove (multsof p) nos) ps...] -- pattern matching
again with the two interpretations/implementations of multsof/remove possible,
one for the trial division, the other - sieve of Eratosthenes. Much more efficiently done,
this time, the both of them.
(we can even add wheels to this, still preserving the ambiguity)
--
for comparison, here's the same functions in KRC-like (well, Haskell) strict pseudocode:
isPrime n lop = null [() | p <- lop, rem n p == 0]
-- OP:
findprimes 2 = [2]
findprimes n = lop ++ [n | isPrime n lop]
where lop = findprimes (n-1)
-- 3rd:
findprimes n | n < 5 = [2,3]
findprimes n = lop ++ [n | n <- [q+1..n], isPrime n lop]
where lop = findprimes q ;
q = floor $ sqrt $ fromIntegral n
-- 2nd:
findprimes n = g 5 9 [2] [3] []
where
g k q a b c
| k > n = a ++ b ++ reverse c
| k == q, [h] <- b = g k q a (h:reverse c) []
| k == q, (h:p:ps) <- b = g (k+2) (p*p) (a++[h]) (p:ps) c
| isPrime k a = g (k+2) q a b (k:c)
| otherwise = g (k+2) q a b c
impl this with generators.... gen's will be able to change their data/definitions inside them... with lexical scope we can set! DECLARED-IN-ADVANCE lexical variables//slots-of-meaning but we can't ADD NEW SLOTS at run time, we must envisage them all in advance...
--
isprime n lop = null [() | p<-lop,rem n p==0] findprimes 2 = [2] findprimes n = lop++[n | isprime n lop] where lop=findprimes(n-1).
measuring (time (length (find-primes 100000))) with your exact code (under #lang racket)
runs for 7 (seven!) seconds in my DrRacket. Not 86. Empirically, it's near quadratic.
--
speaking of painters
and complexities,
using append destroys the algorithmic gains of using sqrt.
your code is just as quadratic as OP's; runs at constant ~4x factor faster (give or take).
set-mcdr! seems unavoidable. or arrays.
-- ("arrays", i.e. dynamic arrays, i.e. list of vectors a-la HAT with fixed size leaf vector blocks, adding more "blocks" as needed, automatically - as we only add elements, never remove, here. or in the general case, arrays with fill-pointer and automatic geometric expansions in the background for O(1) amortized costs... is there an SRFI like that, I wonder?)
---
mcdr--primes.rkt:
#lang racket
(define (is-prime? n lop) ; by Oscar Lopez originally
(define sqrtn (sqrt n))
(let loop ([lop lop])
(cond [(or (empty? lop) (< sqrtn (mcar lop))) true] ; mcar mcons etc. by me
[(zero? (modulo n (mcar lop))) false]
[else (loop (mcdr lop))])))
(define (find-primes n)
(let* ([a (mcons 3 '())]
[b (mcons 2 a)])
(let loop ([p a] [i 5] [d 2] ; d = diff +2 +4 +2 ... ; i +=2 +=4 by me
[c 2]) ; c = count of primes found ; ABSOLUTELY NO
(cond [(> i n) c] ; improvement in speed ?!???
[(is-prime? i (mcdr a))
(set-mcdr! p (mcons i '())) ; (although I measured w/
(loop (mcdr p) (+ i d) (- 6 d) (+ c 1))] ; debugging 'ON'...)
[else (loop p (+ i d) (- 6 d) c )]))))
; (time (length (find-primes 10000))) ; with cons times are in milliseconds
; 10K 156 ; 20K 437 ; 40K 1607 ; 80K 5241 ; 100K 7753 .... n^1.8-1.9-1.7 OP's
; 10K 62 ; 20K 109 ; 40K 421 ; 80K 1217 ; 100K 2293 .... n^1.8-1.9 Oscar Lopez's
; mcons:
(time (find-primes 100000))
; 100K 47 ; 200K 172 ; 1M 1186 ; 2M 2839 ; 3M 4851 ; 4M 7036 .... n^1.25-1.32 this, by me
; 78498 148933 216816
--------
;; Natural Natural -> Boolean
;; Helper function to avoid writing the handful (= 0 (modulo na nb))
(define (divisible na nb) (= 0 (modulo na nb)))
;; Natural ListOfNatural -> Boolean
;; n is the number to check, lop is ALL the prime numbers less than n
(define (is-prime? n lop)
(let ([r (sqrt n)])
(cond [(null? lop) #t]
[(> (car lop) r) #t] ; the only tweak
[(divisible n (car lop)) #f] ; inline or no, no difference
[ else (is-prime? n (cdr lop))])))
;; Natural -> ListOfNatural
(define (find-primes-0 n)
(if (= n 2)
(list 2)
(let ([LOP (find-primes-0 (sub1 n))])
(if (is-prime? n LOP)
(append LOP (list n)) ; traversing the whole list for each new prime
LOP))))
(define (find-primes n)
(let loop ([k 3] ; this candidate
[LOP (list 2)]) ; lop so far
(if (> k n)
LOP
(loop (+ 1 k)
(if (is-prime? k LOP)
(append LOP (list n))
LOP)))))
(time (begin
(display (length (find-primes 100000))) ; 2.10 -> 0.77 ( improvement due to sqrt )
(newline))) ; ( n -> n-2 : 0.70s, 2.25s )
(time (begin
(display (length (find-primes 200000))) ; 6.65 -> 2.35 ( sqrt or no, it's about
(newline))) ; the same ~ n^1.65 )
; code by Oscar Lopez: 0.55, 1.84 ~ n^1.75
my further two variants, running at ~ n^1.3 each:
;; find primes up to and including n, n > 2
(define (find-primes n)
(let loop ( [k 5] [q 9] ; next candidate; square of (car LOP2)
[LOP1 (list 2)] ; primes to test by
[LOP2 (list 3)] ; more primes
[LOP3 (list )] ) ; even more primes, in reverse
(cond [ (> k n)
(append LOP1 LOP2 (reverse LOP3)) ]
[ (= k q)
(if (null? (cdr LOP2))
(loop k q LOP1 (append LOP2 (reverse LOP3)) (list))
(loop (+ k 2)
(* (cadr LOP2) (cadr LOP2)) ; next prime's square
(append LOP1 (list (car LOP2)))
(cdr LOP2) LOP3 )) ]
[ (is-prime? k (cdr LOP1))
(loop (+ k 2) q LOP1 LOP2 (cons k LOP3)) ]
[ else
(loop (+ k 2) q LOP1 LOP2 LOP3 ) ])))
;; n is the number to check, lop is list of prime numbers to check it by
(define (is-prime? n lop)
(cond [ (null? lop) #t ]
[ (divisible n (car lop)) #f ]
[ else (is-prime? n (cdr lop)) ]))
(define (find-primes n)
(cond
((<= n 4) (list 2 3))
(else
(let* ([LOP (find-primes (inexact->exact (floor (sqrt n))))]
[lp (last LOP)])
(local ([define (primes k ps)
(if (<= k lp)
(append LOP ps)
(primes (- k 2) (if (is-prime? k LOP)
(cons k ps)
ps)))])
(primes (if (> (modulo n 2) 0) n (- n 1)) '()))))))
--
__________________________________________________________________________________
isPrime n lop = null lop || not (rem n (head lop) == 0) && isPrime n (tail lop)
= foldr (\p r -> not (rem n p == 0) && r) True lop
= foldr (\b r -> not b && r) True [rem n p == 0 | p <- lop]
= and [ rem n p > 0 | p <- lop]
= all ((> 0) . rem n) lop
= null [() | p <- lop, rem n p == 0]
-- OP:
findprimes 2 = [2]
findprimes n = lop ++ [n | isPrime n lop]
= lop ++ [n | n <- [q+1..n], isPrime n lop]
where lop = findprimes q ; q = n-1
-- 3rd:
findprimes n | n < 5 = [2,3]
findprimes n = lop ++ [n | n <- [q+1..n], isPrime n lop]
where lop = findprimes q ;
q = floor $ sqrt $ fromIntegral n
-- 2nd:
findprimes n = g 5 9 [2] [3] []
where
g k q a b c
| k > n = a ++ b ++ reverse c
| k == q, [h] <- b = g k q a (h:reverse c) [] ---- good
| k == q, (h:p:ps) <- b = g (k+2) (p*p) (a++[h]) (p:ps) c
| isPrime k a = g (k+2) q a b (k:c)
| otherwise = g (k+2) q a b c
__________________________________________________________________________________
findprimes n = g 5 9 [2] [3] []
where
g k q a b c
| k > n = a ++ b ++ reverse c
| k == q, (h:p:ps) <- b = g (k+2) (p*p) (a++[h]) (p:ps) c
| k == q, [h] <- b,
(p:ps) <- reverse c = g (k+2) (p*p) (a++[h]) (p:ps) [] ---- too imperative,
| isPrime k a = g (k+2) q a b (k:c) ---- disparate
| otherwise = g (k+2) q a b c
findprimes n = g 5 9 [2] [3] [] ---- a total mess
where
g k q a b c
| k > n = a ++ b ++ reverse c
| k == q = g (k+2) (p*p) (a++[h]) (p:ps) c2
| otherwise = g (k+2) q a b ([k | isPrime k a] ++ c)
where
(h:t) = b
(p:ps,c2) | null t = (reverse c,[])
| otherwise = (t,c)
