The First N Primes
October 14, 2011
The Sieve of Eratosthenes computes all primes less than some given n. But sometimes you need the opposite, the first n primes, not knowing the largest prime in the set. One method is to calculate the nth prime, using a prime-counting formula, then use the Sieve of Eratosthenes to generate the list. Or instead of calculating the nth prime directly, you can over-estimate the nth prime using the formula Pn < n log (n log n), then trim the list as needed; that’s the natural logarithm to the base e.
Another method is described by Melissa O’Neill in her article The Genuine Sieve of Eratosthenes. Her method is based on sifting each prime, and makes exactly the same sequence of calculations as the regular Sieve of Eratosthenes, but instead of marking bits in a bitarray it uses a priority queue to keep track of composites, one item in the priority queue for each prime; thus, instead of marking all the composite multiples of each sifting prime all at once, in the bitarray, O’Neill’s method keeps a priority queue of the next item in each sifting sequence, inserts a new sifting sequence in the priority queue as each new prime is encountered, and resets the next item in each sifting sequence as it is reached.
Thus, ignoring 2 and the even numbers, O’Neill’s method begins with i = 3 and an empty priority queue. Since 3 is not in the priority queue, it is prime, so it is output, the pair i×i/i×2 = 9/6 is added to the queue, and i is advanced to i + 2 = 5. Now since 5 is not in the queue, it is prime, so it is output, the pair 25/10 is added to the queue, and i is advanced to 7. Now since 7 is not in the queue, it is prime, so it is output, the pair 49/14 is added to the queue, and i is advanced to 9; at this point the queue is [9/6 25/10 49/7]. Now i = 9 is in the queue, so it is composite, and the pair 9/6 is removed from the priority queue and replaced by the pair 15/6, where 15 is calculated as the current composite 9 plus the skip value 6. And so on, until the desired number of primes has been generated. Note that in some cases a composite will be a member of multiple sifting sequences, and each must be replaced; for instance, when i reaches 45, it is a member of the sifting sequences for both 3 and 5, so the priority queue will have pairs 45/6 and 45/10, which should be replaced by 51/6 and 55/10, respectively.
Your task is to implement O’Neill’s algorithm to compute the first n primes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
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]