Postponed Sieve algorithm with start logic

Question

Based on this Python answer by Will Ness, I've been using a JavaScript adaptation for the postponed sieve algorithm from that answer:

function * primes() {
    yield 2;
    yield 3;
    yield 5;
    yield 7;
    const sieve = new Map();
    const ps = primes();
    ps.next() && ps.next();
    for (let p = 3, i = 9; true; i += 2) {
        let s = sieve.get(i);
        if (s !== undefined) {
            sieve.delete(i);
        } else if (i < p * p) {
            yield i;
            continue;
        } else {
            s = 2 * p;
            p = ps.next().value;
        }
        let k = i + s;
        while (sieve.has(k)) k += s;
        sieve.set(k, s);
    }
}

But now that I need to add start point to it, I am struggling to wrap my head around it, for the logic here isn't straightforward.

When start is a prime, I need it to be the first value. And when start is not a prime, I need the sequence to start with the first prime after start.

Will Ness suggested in one of the comments:

You would have to come up with the valid sieve dictionary for the start point. It's easy to get the needed primes - just run this algo up to sqrt(start), then you need to find each core prime's multiple just above (or is it just below?) the start, while minding the duplicates.

Wielding this into reality however is not that simple (or at least for me) :|

Can anybody help with updating this algorithm for such *primes(start) implementation (preferably in JavaScript as above)?

function * primes(start) {

  // how to amend it to support 'start' logic???

}

Conclusion

Following the excellent answer by Will Ness, I decided to share the final code I used, through a public library - primes-generator. All main algorithms can be found in src/sieve.ts.

Jilt

why no comments in the code? is it that trivial and clear what;s going on at each step? what is i? s? p? why no comments in the code? – Cephalopod

@WillNess Because I'm not the author of this adaptation, I copied the JavaScript version from an online resource somewhere in the past. At least, it resembles closely what you wrote, to which I provided link, and that has good documentation ;) – Jilt

interesting, thanks for the clarification. :) I'm actually writing an answer... – Cephalopod

do they give any attribution at all, at that site? – Cephalopod

@WillNess I can no longer retrace where I copied it from, but I remember it was copied exactly as is, without comments. Thank you for looking into it. I was going to put 500 bounty on it just as possible; I think this question well deserves it. – Jilt

I've no objection to the bounty. :) It's actually not that hard though. :) – Cephalopod

thanks for the bounty, and for the chance to write this answer. :) – Cephalopod

Answer 1

(update: added working JS code at the bottom of this answer).

This is the sieve of Eratosthenes:

primes = {2,3,...} \ ⋃(p ← primes) {p², p²+p, ...} = {2} ∪ oddPrimes , oddPrimes = {3,5,...} \ ⋃(p ← oddPrimes) {p², p²+2p, ...}

where \ is set difference (read "minus"), ∪ set union, and ⋃ big union of sets.

An example might illustrate:

{2,3,4,5,6,7,8,9,10,11,12,13,14,15,...}
 \             |
   { 4,  6,  8,| 10,   12,   14,   16,   18, ...}
    \          .
             { 9,      12,      15,      18,    21,  ...}
       \ 
                                               { 25, 30, 35, ...}
          \                                     { 49, 56, 63, ...}
            \                                     { 121, 132, 143, ...}
              \  ........

or for the odd primes,

{3,5,7,9,11,13,15,17,19,21,23,25,27,29,31, ...}
 \                    |
     { 9,      15,    | 21,      27,      33, ...}
   \                          .
                            { 25,            35, ...}
     \                                        { 49, 63, ...}
      \                                         { 121, 143, ...}
       \  ........

The code you refer to in the question implements this approach. At any point in time, when considering a certain candidate, sieve exists in a certain state, as are the rest of the variables in the loop. So we just need to re-create this state directly.

Let's say we're considering i = 19 as the current candidate. At that point we have sieve = { (21, 6) } and p = 5.

This means for a candidate i, sieve contains multiples of all the primes q such that q^2 < i, and p is the next prime after the qs.

Each of the multiples is the smallest one not smaller than i, and there are no duplicates in the sieve. Then it's in a consistent state and can be continued from that point on.

Thus, in pseudocode:

primes() = ..... // as before

primesFrom(start) =
  let
  { primes.next()
  ; ps = takeWhile( (p => p^2 < start) , primes() )
  ; p = primes.next_value()
  ; mults = map( (p => let 
                       { s = 2*p 
                       ; n = (start-p)/s  // integer division NB!
                       ; r = (start-p)%s
                       ; m = p + (r>0 ? n+1 : n)*s
                       }
                       ( m, s) )
                , ps)
  }
  for each (m,s) in mults
    if m in sieve
       increment m by s until m is not in sieve
    add (m,s) to sieve

and then do the same loop as before.

---

As requested, the JS code:

function *primesFrom(start) {
    if (start <= 2) { yield 2; }
    if (start <= 3) { yield 3; }
    if (start <= 5) { yield 5; }
    if (start <= 7) { yield 7; }
    const sieve = new Map();
    const ps = primesFrom(2);
    ps.next();                 // skip the 2
    let p = ps.next().value;   // p==3
    let psqr = p * p;          // p^2, 9
    let c = psqr;              // first candidate, 9
    let s = 6;                 // step value
    let m = 9;                 // multiple

    while( psqr < start )      // must adjust initial state
    {
         s = 2 * p;            
         m = p + s * Math. ceil( (start-p)/s );  // multiple of p
         while (sieve.has(m)) m += s;
         sieve.set(m, s);
         p = ps.next().value;
         psqr = p * p;
    }
    if ( start > c) { c = start; }
    if ( c%2 === 0 ) { c += 1; }
    
    for ( ;  true ;  c += 2 )     // main loop
    {
        s = sieve.get(c);
        if (s !== undefined) {
            sieve.delete(c);      // c composite
        } else if (c < psqr) {
            yield c;              // c prime
            continue;
        } else {                  // c == p^2
            s = 2 * p;
            p = ps.next().value;
            psqr = p * p;
        }
        m = c + s;
        while (sieve.has(m)) m += s;
        sieve.set(m, s);
    }
}

Correctly produces the first 10 primes above 500,000,000 in no time at all on ideone:

Success #stdin #stdout 0.03s 17484KB
500000003,500000009,500000041,500000057,500000069,
500000071,500000077,500000089,500000093,500000099

Apparently it does so with the whopping recursion depth of 5 (five) invocations.

The power of the repeated squaring! or its inverse, the log log operation:

log2( log2( 500000000 )) == 4.85

Cephalopod

Is there any chance you would complete the answer with the updated JavaScript code I included in the question? – Jilt

I need something simple where I just click run and see how it goes. actually, I can do this in ideone.com myself, I see it has JS. :) if I'll have any questions, I'll ask you here. :) --- as for something ready, here's something. – Cephalopod

Can you access this code? – Jilt

Looking forward to your full solution. I have tried following that link to the Haskell-based implementations, even tried to compile it into JavaScript, but without much success so far. I think a lot of people would love seeing it in JavaScript these days ;) – Jilt

I've come to think that sieve algorithm is not as efficient as it is praised. I was able to create a highly-optimized isPrime version, which gives me primes significantly faster than sieve can. – Jilt

Isn't that what I gave you in my comment? It has isPrime implementation + generator. The complete set of tests I have locally - I'm not sure how to include all that. But the idea based on that is trivial. You just generate N numbers and measure time. I tested for generating 10mln primes, got about 5s from sieve, 600ms from isPrime. – Jilt

aha, OK, missed that. could you please modify it to return the first 500 primes after 500,000,000, as my code is doing? does it take 0.04 seconds too? – Cephalopod

I think you've edited your comment.... if it's 10x faster for the 10mln primes, it's great, but it'll become less and less so as we go along. if you test divisibility, the complexity is around n^1.5, and for sieve it's around n^1.0..1.1. so the sieve should win out, eventually. – Cephalopod

So I added the test you asked - generating 50 primes after 50,000,000, and I'm getting about 0.004, i.e. 10 times faster than sieve. – Jilt

for instance, I've modified my code to return first 500 primes after 500 billion, and it took 0.3s on ideone. and for 250 billion, 0.15s. so it grows linearly. I suspect your code won't fair as well. – Cephalopod

See above my test I added as you asked - it goes through those high numbers like they are nothing, outperforming sieve 10 times over. Let me know what you think ;) – Jilt

Generating 50 primes after 500bln takes about 0.7s. I think this is still damn impressive :) – Jilt

so it's now 2x slower than the sieve. – Cephalopod

Yes, but up to about 1bln it is about 10 times faster, and that in itself is interesting. – Jilt

The bounty is now active, so whenever you can finish that sieve update with start logic ;) Cheers! – Jilt

I saw it, but thought you were still working on it, because it makes use of undefined m variable, so the code cannot even run. – Jilt

could you please remove the two preceding comments so as not to confuse the readers? the code runs, and I've edited the code to add the let m = ... line. – Cephalopod

Why did you set m = 9? I do not see 9 having any significance, the initial variable value is not used at all. – Jilt

it doesn't matter. you said it must be defined, so I put some value to it, with let. I don't really know what let does. when m is actually used later on in the code, it is already set to a correct value. – Cephalopod

That's fine, I will be refactoring this code anyhow. Maybe even will publish a library using it, so it is easier for others. I will award the bounty when it becomes possible (22 hours away, for now). Thank you for your help! – Jilt

See the Conclusion note I added to my question. – Jilt

@Jilt great, noted! :) – Cephalopod

@WillNess for large values of start (>2^54) you may experience errors from a loss of precision in the line using "m = p + s * Math.ceil( (start-p)/s )". I think this stems from an intermediary float value only having 53 bits of precision. In a python variant inspired by your code I had to exchange "m = p + s * math.ceil((start - p) / s)" for this conditional modulo operation "m = (start + s - ((start - p) % s)) if ((start - p) % s) else start" to avoid accidentally yielding composites due to incorrect setup. – Bleed

Interesting, thanks for your input. trying this now at ideone, it looks like it would take about a 100 seconds to run there, for 10^18. I wonder how much time did it take on your setup? Were you able to parallelize it perhaps somehow? Thanks. @Bleed – Cephalopod

No parallelism was used! With a start of ~2^54 I was seeing "time to first prime yielded" of about 35 seconds, and with ~2^55 about 50 seconds. I've been using the your prime generator (and a simplistic MillerRabin & Fermat one) to teach myself how to use the unittest module and caught some interesting failures at high start values. @WillNess – Bleed

Answer 2

I have modified Will Ness's amazing answer

1. To allow any start value for printing. As explained below, I think there is no way to avoid it having to compute all the earlier primes, if you want an infinite sequence of higher primes.

2. Adjusted variable names to assist understanding the algorithm.

3. Changed what is stored in the map from 2 times a prime factor, to just the prime factor, to make it easier for readers to follow the algorithm.

4. Moved one part of the code into a subfunction, again for ease of reader understanding.

5. Changed the control flow of the 3-way choice in the middle of the algorithm, and added comments, that simplify understanding. It is probably very slightly slower, because it no longer tests the commonest-true condition first, but it is easier for readers.

function* primeGenerator() {
  yield 2;
  yield 3;
  yield 5;
  yield 7;
  const multiplesWithAFactor = new Map();

  // We only need to consider potential factors  that are themselves prime
  // Start with 3 and get further prime factors on demand from a child instance of ourself.
  let biggestFactorConsidered = 3;
  const childStreamOfPrimes = primeGenerator();
  childStreamOfPrimes.next(); // Discard the 2
  childStreamOfPrimes.next(); // Discard the 3

  let candidate = 7; // We have already yielded 7 as a prime above, so on the first iteration of the while loop we will be testing 9.
  while (true) {
    candidate += 2;

    // Assess candidate, into one of three outcomes: square, non-square multiple, prime. This order is arranged for ease of understanding, but for maximum speed efficiency you should test in the order in Will Ness's version.
    const factorOfCandidate = multiplesWithAFactor.get(candidate);
    
    if (candidate === biggestFactorConsidered * biggestFactorConsidered) {
    // A square, so from now on we will need to consider the next factor, too.
      markNextUnmarkedMultiple(candidate, biggestFactorConsidered);
      biggestFactorConsidered = childStreamOfPrimes.next().value;
    } else if (factorOfCandidate) {
    // A non-square multiple. Can forget that fact now, and instead mark the next unmarked multiple of the same prime factor
      multiplesWithAFactor.delete(candidate);
      markNextUnmarkedMultiple(candidate, factorOfCandidate);
    } else {
    // Prime
      yield candidate;
    }
  }

  // This is a subfunction for ease of understanding, but for maximum efficiency you should put this in the while loop above, and have the "yield candidate" avoid calling it, via a "continue" statement.
  function markNextUnmarkedMultiple(candidate, factor) {
    let nextMultiple = candidate + 2 * factor;
    while (multiplesWithAFactor.has(nextMultiple)) {
      nextMultiple += 2 * factor;
    }
    multiplesWithAFactor.set(nextMultiple, factor);
  }
}

const maxItems = 20;
const minimumPrintablePrime = 5e8;
console.log("Running");
const primeObject = primeGenerator();
let count = 0;
do {
  const prime = primeObject.next().value;
  if (prime > minimumPrintablePrime) { // If you don't like seeing this process of reading non-printed primes out here in the parent code, you can do it within the primeGenerator itself. Either way, _someone_ has to call the prime generator for all primes up to the square root of the starting value.
    console.log(prime);
    count++;
  }
} while (count < maxItems);

The algorithm is astonishingly economical on storage

It uses several instances of the prime generator. You call the "main" one, which runs along examining candidates, while storing a set of future candidates that it knows are multiples of smaller prime factors.

At any one time, when considering candidate N, its stored map consists of an entry for each prime factor up to sqrt(N) but stored in the map at the next multiple of the prime factor which the algorithm will reach, as it explores forward from the current candidate.

Whenever the candidate is present in this map of multiples, it can reject that candidate. The map tells us, for that multiple, what was the prime factor which told us that number was a multiple. For example, 15 is marked as a multiple of 3. When the algorithm comes to 15, it realises it is a multiple of 3, so it rejects the 15.

From now on, it no longer needs to remember that 15 is a multiple, because 15 will never be a candidate again in this instance of the algorithm. It can delete the 15 entry from its map of multiples, and replace it with the next multiple of 3. However the very next multiple will always be even, because all candidates (e.g. 15) and all factors (e.g. 3) are odd. Therefore it only needs to tell the map that 15 + 2x3 is a multiple, i.e. it can step forward by 2 x the factor. Moreover, if the resulting number, 21, is already marked as a multiple, it need not mark that number: it can advance by a further 2x3 to 15 + 2x3 + 2x3, etc, until it finds a number that is not already marked as a multiple.

Cleverly, this process ensures that every factor (e.g. 3) remains marking a multiple in the map forever. To assess a candidate N, the map will only contain one entry for every prime up to sqrt(N).

When the candidate rises above the square of the largest factor considered so far, the algorithm needs to get the next factor. It simply uses a second instance of itself to get the next prime in the sequence. I was at first worried about this creating vast memory demands, but that is not so. Assessing a number ~2^64, would need all the second instance for ~2^32, which in turn would call a third instance for ~2^16, and so on. Only a handful of instances of the generator need to be created even for primes of enormous size, and if an instance's map has factors up to F, the next instance only needs factors up to sqrt(F), which rapidly becomes small.

Even the Ness algorithm still needs to build a map of all the factors up to sqrt(N)

It needs the map to contain all those factors, so that it can correctly reject candidates around N.

But it also will _eventually_ need the factors between sqrt(N) and N

Because only that way can it be confident of numbers from N to N^2 being prime.

My conclusion

I think it needs to iterate through, in order to work. I can't see a way of starting it at an arbitrary candidate, e.g. 10^8, without pre-filling the map with (all the) primes up to 10^4.

The most efficient way of doing that pre-filling, is by running the algorithm itself, which is in fact what it does if you simply ask it to yield (but do not print) all the primes up to 10^8.

My suspicion is that the algorithm is so efficient, that the best way to get it into position to start at an arbitrary candidate, e.g. 1e8, is by running the algorithm itself, i.e. there is no shortcut. Remarkable situation!

This seems to be confirmed in Will's updated answer, in which the algorithm calls itself to prefill the sieve. At the end of the prefilling, the sieve contains one entry for each prime up to sqrt(start). My version is much slower, but this is not because it collects more primes than Will's; it is because it (a) tests for the 3 conditions in a different order (to make it easier for humans to follow) and (b) extracts some code into a subfunction (again to make it easier to understand). Obviously Will's approach is better for production code!

Lexicographer

@will Ness - thank you! Is my revised answer correct? You seem to be very expert in Prime generation. Am I right that, if you want a generator that goes on forever, it needs to begin with knowing, explicitly or implicitly, all primes up to the start value? – Lexicographer

implicitly every true thing in the Universe is True. just by virtue of it being true. but we're in the business of the explicit knowledge. the sieve of Eratosthenes is primes = [2,3,...] \ U [[p*p,p*p+p,...] for p in primes]. or more verbosely primes = nats2 \ U composites, composites = map multiples primes, multiples(p) = [p*p,p*p+p,...] -- the multiples of each prime p start being enumerated from p*p. so to find a prime q only the primes up to sqrt(q) need to be known explicitly. the knowledge of all the primes up to q is implicit in the primes up to sqrt(q). – Cephalopod

so all we really need to know is that the first prime is 2. everything else follows. @Eureka – Cephalopod

I've posted an answer here as well. in particular it does describe a shortcut, and it shows that there's no consideration of any factors in the sieve. that would be the way of a trial division sieve, not the sieve of Eratosthenes. :) – Cephalopod

I give full attribution in the Python answer linked in the question. in particular, the whole dict thing isn't mine at all. my contribution was the recursive call to primes(). :) the length of that nested "tower" of self-invocations is O(log log n), because of the repeated taking of the square root. (which was not a red herring after all, huh. :) :) ) – Cephalopod

Oh no, you just let it run until it hits the starting threshold. this isn't what was asked for in the question, nor described in the answer. :) I'll still have to finish that JS code. dang. :) I'll edit your answer to add a starting notice. it shows 4 invocations, and full 5 seconds gap before it gets to that starting point. the whole point is that it should take near 0 time before it starts. – Cephalopod

"I can't see a way of starting it at an arbitrary candidate, e.g. 10^8, without pre-filling the map with (all the) primes up to 10^4" correct. "The most efficient way of doing that pre-filling, is by running the algorithm itself," incorrect. consider a simpler question: what is the smallest odd number above 1E8? what you say is the same as saying, let's compute {x=1; do {x+=2} while (x<1E8)}. which is correct, but then you say it's the most efficient way. it's not. :) – Cephalopod

FYI I've added a working JS code to my answer. there is a shortcut after all. :) – Cephalopod

It's brilliant Will! But the line "p = ps.next().value;" you are indeed calling the algorithm itself, to pre-fill "sieve" with a value for each prime up to sqrt(N). If you don't believe me, console.log sieve after that sieve-initialization phase. It has an entry for each prime up to sqrt(N), and it got it by calling the child instance of the same function. – Lexicographer

so why "but"? what's new about it? but your answer says it must compute all previous primes. this is just plain wrong. much of the rest of the stuff in there is wrong a well. question: what's the smallest odd multiple of 7 above 10^8? my answer: ceil((10^8-7)/14)*14+7. your answer is: for (i=7; i < 10^8; i+=14) {}; print(i);. which one is waaaay slower, do you think? which one is a shortcut to find the answer instantaneously? there is a shortcut, and my answer does not "confirm" it doesn't exist, it contradicts that claim, and actually shows the shortcut. – Cephalopod