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.

About these ads

Pages: 1 2

5 Responses to “The First N Primes”

  1. Mike said

    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))
    
    
  2. CyberSpace17 said

    My C++ solution

    http://ideone.com/CWsQx

  3. CyberSpace17 said

    http://ideone.com/B0Bpy

    my solution again with an amendment to line 18

  4. David said

    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   509 
    
  5. David said

    This 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]
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 647 other followers

%d bloggers like this: