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 articleThe 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.

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…

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):

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:

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:

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:

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.

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:

@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