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.
Pages: 1 2
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]