The First N Primes
October 14, 2011
Our solution pre-computes the first element of the sieve, thus avoiding a test for an empty priority queue at each step through the algorithm:
(define (n-primes n)
(let ((pq (pq-insert lt? (cons 9 6) pq-empty)))
(let loop1 ((i 5) (pq pq) (ps (list 3 2)) (k 2))
(cond ((= n k) (reverse ps))
((< i (car (pq-first pq)))
(let* ((c (* i i)) (s (+ i i))
(pq (pq-insert lt? (cons c s) pq)))
(loop1 (+ i 2) pq (cons i ps) (+ k 1))))
(else (let loop2 ((pq pq))
(if (< i (car (pq-first pq)))
(loop1 (+ i 2) pq ps k)
(let* ((c (car (pq-first pq)))
(s (cdr (pq-first pq)))
(pq (pq-rest lt? pq)))
(loop2 (pq-insert lt? (cons (+ c s) s) pq))))))))))
The priority queue is ordered by increasing i×i, and then by increasing i×2 in case of a tie:
(define (lt? a b)
(or (< (car a) (car b))
(and (= (car a) (car b))
(< (cdr a) (cdr b)))))
Here is a list of the first 168 primes:
> (n-primes 168)
(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 547 557 563 569 571 577 587 593 599 601 607 613
617 619 631 641 643 647 653 659 661 673 677 683 691 701 709
719 727 733 739 743 751 757 761 769 773 787 797 809 811 821
823 827 829 839 853 857 859 863 877 881 883 887 907 911 919
929 937 941 947 953 967 971 977 983 991 997)
We used priority queues from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/61r9U0eR.
Programming Praxis 
Python 3 version. Uses heapq from the standard library for the priority queue.
I began odd_number at 5 to eliminate a test for when ‘next_composite’ isn’t set yet
from heapq import * from itertools import count, islice def oneill(): yield 2 yield 3 queue = [] next_composite, skip = 9, 6 for odd_number in count(5, 2): if odd_number < next_composite: yield odd_number heappush(queue, (odd_number**2, 2*odd_number)) else: while odd_number == next_composite: next_composite, skip = heappushpop(queue, (next_composite+skip, skip)) def first_n_primes(n): return list(islice(oneill(), n))My C++ solution
http://ideone.com/CWsQx
http://ideone.com/B0Bpy
my solution again with an amendment to line 18
Factor version, uses Priority Queue library and lazy lists.
USING: kernel math heaps sequences arrays lists lists.lazy locals io formatting ; IN: oneil-sieve : o'neil-init ( -- heap ) 6 9 <min-heap> [ heap-push ] keep ; : heap-top ( pq -- n ) heap-peek nip ; : isprime? ( pq n -- ? ) [ heap-top ] dip = not ; : push-prime ( pq n -- ) [ 2 * ] [ sq ] bi rot heap-push ; :: advance ( pq n -- ) [ pq heap-top n = ] [ pq heap-pop over + pq heap-push ] while ; : o'neil-step ( pair -- pair ) first2 2 + [ 2dup isprime? not ] [ [ advance ] 2keep 2 + ] while [ push-prime ] 2keep 2array ; : o'neil-sieve ( -- lazy-list ) o'neil-init 3 2array [ o'neil-step ] lfrom-by [ second ] lazy-map ; : all-primes ( -- lazy-list ) o'neil-sieve 2 swons ; CONSTANT: items/line 10 : nl-cond ( counter list -- counter list ) [ 1 + ] dip over items/line = [ [ drop 0 ] dip nl ] when ; : .primes ( lazy-list -- ) 0 swap [ dup nil? ] [ unswons "%5d " printf nl-cond ] until drop 0 = not [ nl ] when ;Test:
( scratchpad ) all-primes [ 512 < ] lwhile .primes 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 509This method is better suited to Prolog than the standard sieve, due to the fact that arrays do not exist in that language, and generators are pretty much how things are done, so without further ado, in Prolog…
% O'neil sieve ?- use_module(library(heaps)). prime(2). prime(N) :- prime_heap(N, _). prime_heap(3, H) :- list_to_heap([9-6], H). prime_heap(N, H) :- prime_heap(M, H0), N0 is M + 2, next_prime(N0, H0, N, H). next_prime(N0, H0, N, H) :- \+ min_of_heap(H0, N0, _), N = N0, Composite is N*N, Skip is N+N, add_to_heap(H0, Composite, Skip, H). next_prime(N0, H0, N, H) :- min_of_heap(H0, N0, _), adjust_heap(H0, N0, H1), N1 is N0 + 2, next_prime(N1, H1, N, H). adjust_heap(H0, N, H) :- min_of_heap(H0, N, _), get_from_heap(H0, N, Skip, H1), Composite is N + Skip, add_to_heap(H1, Composite, Skip, H2), adjust_heap(H2, N, H). adjust_heap(H, N, H) :- \+ min_of_heap(H, N, _). prime_seq(Test, Prime) :- prime(P), (call(Test, P) -> Prime = P ; !, false). upto_1000(X) :- X < 1000. ?- setof(P, prime_seq(upto_1000, P), S), length(S, L). L = 168, S = [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,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701, 709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911, 919,929,937,941,947,953,967,971,977,983,991,997]Several comments or improvements as follows:
1. O’Neill introduced the idea of Wheel Factorization, in her sample case up to the 48 gap value entries for the 2/3/5/7 wheel but acknowledged that bigger wheels could be used, although she underestimated their value. A 2/3/5/7/11/13/17/19 wheel offers a reduction in execution time to about 58% as compared to using her wheel for the primes up to 100 million.
2. O’Neill didn’t realize it in her referenced paper but in later work realized that the work can be reduced by postponing the adding of a prime composite series to the map/PriorityQueue/HashTable until the range has increased to the square of the current base prime value. This makes the versions in her paper using a (Binary Tree) Map or Priority Queue faster since they add a (log n) term to the performance which would change to the log of the square root of n or a constant factor of one half of log n. It doesn’t help a HashTable implementation as much with their amortized O(1) performance, but it greatly reduces the memory requirements for both from O(n / log n) to O(square root of n / log n).
3. The Priority Queue as used is fine for smaller prime ranges but it does add this log n to the performance making it O(n log n log log n) from O(n log log n); use of a HashTable with amortized O(1) performance gives it the original performance at the cost of the constant factor average overhead of handling the hashing. Even if this constant factor was 10, the HashTable would be better for sieve ranges over about 500 million (the factor may be much less than that depending on the implementations). Most languages offer a standard library function for a HashTable/Dictionary but not Scheme R5RS, although that is correcting in R6RS. Python and JavaScript make it easy because every object is a HashTable (ie. built in to the languages).
4. Incremental sieves based on these are always much slower than incremental sieves based on page segmented arrays due to these constant factor overheads in processing the Map/Priority Queue/HashTable, as sieving arrays consists of nothing but very tight inner culling loops.
My Python version with a huge wheel factorization, parameterized for easy experimentation (wheel larger than this aren’t recommended due to excessive memory use, but one can try smaller to find relative improvements):
FRSTBP = 23; WHLPRMS = [ 2, 3, 5, 7, 11, 13, 17, 19 ] # base wheel primes def makeGAPS(): whlsz = reduce(lambda o, x: o * x, WHLPRMS, 1) # length of the repeating wheel whlhts = reduce(lambda o, x: o * (x - 1), WHLPRMS, 1) # number of candidates in wheel cndptrn = [True] * (whlsz + 1) # wheel size plus one extra for o/f for p in WHLPRMS: strt = p * p - FRSTBP while strt < 0: strt += p cndptrn[strt::p] = [False] * ((len(cndptrn) - strt - 1) // p + 1) # cull for p def getnxtcndndx(i): j = i + 1 while not cndptrn[j]: j += 1 return j def getgaps(): i = 0; for h in xrange(whlhts): j = getnxtcndndx(i) yield j - i; i = j return bytearray(getgaps()) GAPS = makeGAPS() def primes(): def wheel_prime_pairs(): yield (FRSTBP, 0) # avoid initialization race condition bps = wheel_prime_pairs() # additional primes supply, recursion required but not deep p, pi = next(bps); q = p * p # adv to get 11 sqr'd is 121 as next square to put nxtCmpsts = {}; n = FRSTBP + GAPS[0]; ni = 0 # into sieve dict; init cndidate, wheel ndx if len(GAPS) > 1: ni = 1 while True: # Python does not recurse well do not required here if n not in nxtCmpsts: # is not a multiple of previously recorded primes if n < q: yield (n, ni) # n is prime with wheel modulo index else: # n == p * p: put next cull position on wheel nxt = q + p * GAPS[pi]; pi += 1 # advance wheel index if pi >= len(GAPS): pi = 0 while nxt in nxtCmpsts: # ensure each entry is unique by wheel nxt += p * GAPS[pi]; pi += 1 # advance wheel index if pi >= len(GAPS): pi = 0 nxtCmpsts[nxt] = (p, pi) # next non-marked multiple of a prime p, pi = next(bps); q = p * p # advance next prime and prime square to compare else: pc, pci = nxtCmpsts.pop(n) # move current cull position up the wheel nxt = n + pc * GAPS[pci]; pci += 1 # advance wheel index if pci >= len(GAPS): pci = 0 while nxt in nxtCmpsts: # ensure each entry is unique by wheel nxt += pc * GAPS[pci]; pci += 1 # advance wheel index if pci >= len(GAPS): pci = 0 nxtCmpsts[nxt] = (pc, pci) # next non-marked multiple of a prime n += GAPS[ni]; ni += 1 # advance on the wheel if ni >= len(GAPS): ni = 0 for p in WHLPRMS: yield p # root primes for p, pi in wheel_prime_pairs(): yield p # strip out indexes def first_n_primes(n): count = 0; gen = primes() while True: if count >= n: break count += 1; yield gen.next() numPrimes = 10000000 # testing from time import clock elpsd = -clock() lastPrime = 0 for p in first_n_primes(numPrimes): lastPrime = p elpsd += clock() print("The " + str(numPrimes) + "th prime is " + str(lastPrime) + " in " + str(int(elpsd * 1000)) + " milliseconds.")The above code generates the following output (the couple of seconds used to generate the huge wheel isn’t included in the time):
it performs very close to its theoretical performance for a large range of primes.
The following code is an improvement to the above code in that it reduces the overhead time to generate the GAPS byte array (replaces “makeGAPS()”) by very nearly a factor of four from well over two seconds to well under one second on my machine; It won’t affect the timing results as the overhead time is excluded from the timing results:
FRSTBP = 23; WHLPRMS = [ 2, 3, 5, 7, 11, 13, 17, 19 ] # base wheel primes def makeGAPS(): whlsz = reduce(lambda o, x: o * x, WHLPRMS, 1) # length of the repeating wheel whlhts = reduce(lambda o, x: o * (x - 1), WHLPRMS, 1) # number of candidates in wheel cndptrn = [True] * (whlsz + 1) # wheel size plus one extra for o/f for p in WHLPRMS: strt = p * p - FRSTBP while strt < 0: strt += p cndptrn[strt::p] = [False] * ((len(cndptrn) - strt - 1) // p + 1) # cull for p def getnxtcndndx(i): j = i + 1 while not cndptrn[j]: j += 1 return j def getgaps(): i = 0; for h in xrange(whlhts): j = getnxtcndndx(i) yield j - i; i = j return bytearray(getgaps()) GAPS = makeGAPS()Sorry, I selected the wrong version above, and submitted the original version:
The following code is an improvement to the above code in that it reduces the overhead time to generate the GAPS byte array (replaces “makeGAPS()”) by very nearly a factor of four from well over two seconds to well under one second on my machine; It won’t affect the timing results as the overhead time is excluded from the timing results:
FRSTBP = 23; WHLPRMS = [ 2, 3, 5, 7, 11, 13, 17, 19 ] # base wheel primes def makeGAPS(): gaps = bytearray([1]) # initial gap set for p in WHLPRMS: def whlgapsitr(): frstprm = p + gaps[0] n = frstprm; ni = 1; ag = 0; nxtcmpst = p * p; ci = 0 for i in xrange(len(gaps) * p): if ni >= len(gaps): ni = 0 g = gaps[ni]; ni += 1 n += g; ag += g if n < nxtcmpst: yield ag; ag = 0 else: nxtcmpst += p * gaps[ci]; ci += 1 gaps = bytearray(whlgapsitr()) return gaps GAPS = makeGAPS()Since this website likes Scheme code, here is a version of the Richard Bird’s sieve from Epilogue of the same Melissa E. O’Neill artlcle. Although it is quite a bit slower than more imperative algorithms or even incremental sieves using a Priority Queue or HashTable, it is acceptably fast for smaller ranges. It can be sped up by using tree folding (https://www.haskell.org/haskellwiki/Prime_numbers#Tree_merging) and/or gains a constant factor in speed using wheel factorization. The code is as follows:
(define (birdPrimes) (define (mltpls p) (define pm2 (* p 2)) (let nxtmltpl ((cmpst (* p p))) (cons cmpst (lambda () (nxtmltpl (+ cmpst pm2)))))) (define (allmltpls ps) (cons (mltpls (car ps)) (lambda () (allmltpls ((cdr ps)))))) (define (merge xs ys) (let ((x (car xs)) (xt (cdr xs)) (y (car ys)) (yt (cdr ys))) (cond ((< x y) (cons x (lambda () (merge (xt) ys)))) ((> x y) (cons y (lambda () (merge xs (yt))))) (else (cons x (lambda () (merge (xt) (yt)))))))) (define (mrgmltpls mltplss) (cons (car (car mltplss)) (lambda () (merge ((cdr (car mltplss))) (mrgmltpls ((cdr mltplss))))))) (define (minuscmpsts n cmps) (if (< n (car cmps)) (cons n (lambda () (minuscmpsts (+ n 2) cmps))) (minuscmpsts (+ n 2) ((cdr cmps))))) (letrec ((cmpsts (lambda () (mrgmltpls (allmltpls (oddprms))))) (oddprms (lambda () ;; prevent initialization race by presetting 3... (cons 3 (lambda () (minuscmpsts 5 (cmpsts))))))) (cons 2 (lambda () (oddprms))))) (define (nthPrime n) (let nxtprm ((cnt 0) (ps (birdPrimes))) (if (< cnt n) (nxtprm (+ cnt 1) ((cdr ps))) (car ps)))) (nthPrime 1000000)The above code may need a description of how it works other than stating that “it follows the Richard Bird algorithm from the article”:
1. As Scheme doesn’t have a built-in Stream, the above code implements Co Inductive Streams (CIS’s) (http://c2.com/cgi/wiki?CoinductiveDataType) using ‘cons’, ‘car’ and ‘cdr’ from lists, with the difference being that the second element of each pair doesn’t contain a link to the next element of the list but rather a lambda function with no arguments (a “thunk”) that is the computation to obtain the next pair. In a similar means one could use ‘delay’ and ‘force’ to produce a lazy list or lazy stream, but this algorithm does’t need the extra overhead implied by the memoization of ‘force’ in remembering the value so it isn’t computed again on repeated calls, as the above algorithm only passes through each stream once. The difference in its use once constructed is that in order to obtain the next element, instead of calling (cdr element) directly, one calls ((cdr element)) (double brackets) in order to cause the thunk to be applied. As used, these are infinite CIS’s as the algorithm produces an “infinite” stream of primes (or at least until the computations start to overflow – quite quickly for 32-bit versions of Scheme; not for a very long time for 64-bit versions).
2. There are several streams used as follows: an allmltpls stream of streams of all the multiples for each base prime, which is merged into a single stream of composites with repeat composites eliminated (important for the “minus” step), the final output stream of primes, and another separate stream of base primes starting at 3. It is this last requirement that means that ‘cmpsts’ and ‘oddprms’ must be defined recursively (‘letrec’) as they call each other since ‘cmpsts’ uses a separate source of ‘oddprms’ and ‘oddprms’ must use ‘cmpsts’ in the ‘minuscmpsts’ operation.
3. No implementation of a fold right operation is required as that has been wrapped into the ‘mrgmltpls’ operation; no implementation of a mapping operation is required to produce all the multiples of each base prime as that has been combined into the ‘mltpls’ production as a constant addition of twice the prime value series (a ‘scan right’ operation in effect).
4. Note that the primes generator does not produce a list although the output CIS could be passed to a function that could turn CIS’s into conventional Scheme lists; For most purposes that is not necessary as here in showing the last element(s) of the stream up to a limit – direct processing of the stream is all that is necessary.
As per the O’Neil article, this algorithm is somewhat slow as to asymptotic computational complexity at O(n (sqrt n) (log log n) / (log n)2) compared to straight imperative array based at O(n log log n) or even the Map/Priority Queue based versions at O(n (log n) (log log n)), but it is quite interesting, not hard to implement (no external implementations of Priority Queue or HashTable necessary), and will keep up to those implementations up to ranges between about 10 million and 100 million. Tree folding will decrease the run times for a give range as it reduces the depth of the merging to about half with most of the operations in the left (shorter) end so as to extend the useful range of these types of sieve. Adding wheel factorization will reduce the execution time by a constant factor of about 4, but of course can be applied to all of the different sorts of algorithms.
@WillyGood: If you want to use streams, you might want to look at SRFI-41 on the Essays page. I am the author of SRFI-41.
A few further notes:
1. I’ve been doing some work comparing performance of various Scheme implementations and find Gambit-C or Chicken to be about the fastest as well as being fairly fully featured and easy to obtain (I cannot afford or access the actual Chez Scheme compiler – only using the free Petit Chez Scheme byte code interpreter based version – so have no way of comparing its performance), and found some interesting results: Optimized Gambit with “(not safe)” and “(standard-bindings)” declarations “blows Chicken away” for highly recursive continuation style algorithms such as this at about eight times faster than maximally optimized Chicken; however Chicken is over twice as fast implementing very fast tail-call optimized loops such as used in the regular array based sieves – this is likely due to the “Cheney-based” tail cal optimization used by Chicken. Thus Gambit is a better choice for this algorithm and Chicken is a better choice for algorithms such as sieves with very tight recursive loops, where Chickens speed in optimizing these is very close to optimized imperative ‘C’ code; it would be interesting to find out how the full Chez Scheme compiler performs in these.
2. I have learned that internal (define …)’s are actually equivalent to the (letrec …) form used, as internal ”define’s” are mutually recursive (within limits as long as they don’t modify or depend on defined variables just as for ‘letrec”‘s), but internal (define … ) syntax is likely easier to read and understand; I have thus changed the code to:
(define (birdPrimes) (define (mltpls p) (define pm2 (* p 2)) (let nxtmltpl ((cmpst (* p p))) (cons cmpst (lambda () (nxtmltpl (+ cmpst pm2)))))) (define (allmltpls ps) (cons (mltpls (car ps)) (lambda () (allmltpls ((cdr ps)))))) (define (merge xs ys) (let ((x (car xs)) (xt (cdr xs)) (y (car ys)) (yt (cdr ys))) (cond ((< x y) (cons x (lambda () (merge (xt) ys)))) ((> x y) (cons y (lambda () (merge xs (yt))))) (else (cons x (lambda () (merge (xt) (yt)))))))) (define (mrgmltpls mltplss) (cons (car (car mltplss)) (lambda () (merge ((cdr (car mltplss))) (mrgmltpls ((cdr mltplss))))))) (define (minuscmpsts n cmps) (if (< n (car cmps)) (cons n (lambda () (minuscmpsts (+ n 2) cmps))) (minuscmpsts (+ n 2) ((cdr cmps))))) (define (cmpsts) (mrgmltpls (allmltpls (oddprms)))) ;; internal define's are mutually recursive (define (oddprms) (cons 3 (lambda () (minuscmpsts 5 (cmpsts))))) (cons 2 (lambda () (oddprms)))) (define (nthPrime n) (let nxtprm ((cnt 0) (ps (birdPrimes))) (if (< cnt n) (nxtprm (+ cnt 1) ((cdr ps))) (car ps)))) (nthPrime 1000000)Here is the tree folding/merging version in Scheme, showing how easy it is to add the folding with the “pairs” function merged into the “mrgmltpls” function:
(define (treemergePrimes) (define (mltpls p) (define pm2 (fx* p 2)) (let nxtmltpl ((cmpst (fx* p p))) (cons cmpst (lambda () (nxtmltpl (fx+ cmpst pm2)))))) (define (allmltpls ps) (cons (mltpls (car ps)) (lambda () (allmltpls ((cdr ps)))))) (define (merge xs ys) (let ((x (car xs)) (xt (cdr xs)) (y (car ys)) (yt (cdr ys))) (cond ((fx< x y) (cons x (lambda () (merge (xt) ys)))) ((fx> x y) (cons y (lambda () (merge xs (yt))))) (else (cons x (lambda () (merge (xt) (yt)))))))) (define (pairs mltplss) (let ((tl ((cdr mltplss)))) (cons (merge (car mltplss) (car tl)) (lambda () (pairs ((cdr tl))))))) (define (mrgmltpls mltplss) (cons (car (car mltplss)) (lambda () (merge ((cdr (car mltplss))) (mrgmltpls (pairs ((cdr mltplss)))))))) (define (minusstrtat n cmps) (if (fx< n (car cmps)) (cons n (lambda () (minusstrtat (fx+ n 2) cmps))) (minusstrtat (fx+ n 2) ((cdr cmps))))) (define (cmpsts) (mrgmltpls (allmltpls (oddprms)))) ;; internal define's are mutually recursive (define (oddprms) (cons 3 (lambda () (minusstrtat 5 (cmpsts))))) (cons 2 (lambda () (oddprms))))The above code is about three times faster than the Richard Bird version for the same practical range (empirical performance).
@programmingpraxis, your suggestion to use SRFI-41 is all right, but I was already aware of it and decided not to use it, as mentioned in my post, since this algorithm doesn’t require the extra overhead of the memoization implied in using ‘delay’ and ‘force’ to build the “thunk” in the cdr of each pair; also this algorithm didn’t require any of the “higher-order functions” of that SRFI. This is similar to implementing these same algorithms in F# which doesn’t have a built in “lazy list” (one can inport the “PowerPack” for that, if they really need it), but for the same reasons I elected not to do that and used Co-Inductive Streams, just as here
Melissa O’Neill may have proposed the use of a priority queue in 2009, but I first saw the idea as an exercise in Donald Knuth’s The Art of Computer Programming, Vol 3 Exercise 15 of part 5.2.3. The description of the method is in the answers to exercices at pages 641-642 of the 1998 edition.