A simple sieve in Phobos?

[bearophile]

There is a efficient Sieve implementation in C++ here:
http://code.activestate.com/recipes/577966-even-faster-prime-generator/?in=lang-cpp

There are of course far faster implementations, but its 
performance is not bad, while being still simple and quite short.

A D implementation (it's not a final version because the global 
enums should be removed and replaced with run-time variables or 
template arguments. And something different from just counting 
must be added):


import std.stdio, std.algorithm, std.typecons;

enum uint MAX_N = 1_000_000_000U;
enum uint RT_MAX_N = 32_000; // square of max prime under this 
limit should exceed MAX_N.
enum uint B_SIZE = 20_000;   // not sure what optimal value for 
this is;
                             // currently must avoid overflow 
when squared.

// mod 30, (odd) primes have remainders 1,7,11,13,17,19,23,29.
// e.g. start with mark[B_SIZE/30]
// and offset[] = {1, 7, 11, ...}
// then mark[i] corresponds to 30 * (i / 8) + b - 1 + offset[i % 
8].
Tuple!(uint, size_t, uint) calcPrimes() pure nothrow {
    // Assumes p, b odd and p * p won't overflow.
    static void crossOut(in uint p, in uint b, bool[] mark)
    pure nothrow {
        uint si = (p - (b % p)) % p;
        if (si & 1)
            si += p;
        if (p ^^ 2 > b)
            si = si.max(p ^^ 2 - b);

        for (uint i = si / 2; i < B_SIZE / 2; i += p)
            mark[i] = true;
    }

    uint pCount = 1; uint lastP = 2;
    // Do something with first prime (2) here.

    uint[] smallP = [2];

    bool[B_SIZE / 2] mark = false;
    foreach (immutable uint i; 1 .. B_SIZE / 2) {
        if (!mark[i]) {
            pCount++; lastP = 2 * i + 1;
            // Do something with lastP here.

            smallP ~= lastP;
            if (lastP ^^ 2 <= B_SIZE)
                crossOut(2 * i + 1, 1, mark);
        } else
            mark[i] = false;
    }

    for (uint b = 1 + B_SIZE; b < MAX_N; b += B_SIZE) {
        for (uint i = 1; smallP[i] ^^ 2 < b + B_SIZE; ++i)
            crossOut(smallP[i], b, mark);
        foreach (immutable uint i; 0 .. B_SIZE / 2) {
            if (!mark[i]) {
                pCount++; lastP = 2 * i + b;
                // Do something with lastP here.

                if (lastP <= RT_MAX_N)
                    smallP ~= lastP;
            } else
                mark[i] = false;
        }
    }

    return tuple(pCount, smallP.length, lastP);
}

void main() {
    immutable result = calcPrimes;

    writeln("Found ", result[0], " primes in total.");
    writeln("Recorded ", result[1], " small primes, up to ", 
RT_MAX_N);
    writeln("Last prime found was ", result[2]);
}


Is it a good idea to add a simple but reasonably fast Sieve 
implementation to Phobos? I have needed a prime numbers lazy 
range, and a isPrime() function numerous times. (And as for 
std.numeric.fft, people that need even more performance will use 
code outside Phobos).

Bye,
bearophile