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???
}
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.
(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