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.

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 612 other followers

%d bloggers like this: