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

I can't understand whether or not this is a sieve of atkin...


On Tuesday, 18 March 2014 at 14:23:32 UTC, bearophile wrote:
> 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

I've had good luck with Stephan Brumme's block-wise sieve[1], also based on the Sieve of Eratosthenes. It's best when calculating several blocks in parallel, of course. But it's still pretty fast even when used sequentially.

[1]: http://create.stephan-brumme.com/eratosthenes/

On Tuesday, 18 March 2014 at 15:54:23 UTC, Andrea Fontana wrote:
> I can't understand whether or not this is a sieve of atkin...
>
>

The link says 'A very quick (segmented) sieve of Eratosthenes'

On Tuesday, 18 March 2014 at 14:23:32 UTC, bearophile wrote:
> 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).

My new job has finally given me permission to continue some work I was doing for Phobos, which includes an implementation of RSA. RSA uses primes quite heavily, so an isPrime() method is part of that. I think I templated it to accept BigInt or regular ints, but even if not, that should be an easy addition.

I could break out that and a few other basic math functions/algorithms into its own small pull request, if desired?

Chris Williams:

> so an isPrime() method is part of that. I think I templated it to accept BigInt or regular ints, but even if not, that should be an easy addition.
>
> I could break out that and a few other basic math functions/algorithms into its own small pull request, if desired?

I suggest a Phobos module named "combinatorics" (or just "combs"?). It's not meant to be a complete library of combinatorics algorithms, nor to contain the most optimized algorithms around. It's meant to be a collection of efficient but sufficiently short implementations of the few algorithms/ranges you need most often. Everything else, including the most efficient code, I my opinion should be left to specialized numerical libraries external to Phobos (or could be added later if Phobos gains more developers).

I think the most commonly useful functions are:

- A lazy range (a simple unbounded segmented Sieve) that generates primes numbers very quickly, in a given range, or from 2;
- A isPrime() function. Probably it should cache some of its computations.

- A function to compute the GCD on ulongs/longs/bigints is useful.
(Issues 4125 and 7102).

- An efficient and as much as possibly overflow-safe binomial(n,k) that returns a single number.

I'd also like permutations/combinations/pairwise ranges (Phobos already has a permutations, but it's designed on the legacy C++ style of functions, so it's not good enough).
(See ER issue 6788 for pairwise. A the moment I can't find my Bugzilla ER entry for permutations/combinations, but you can see the good API for the permutations/combinations ranges in the code I have written here: http://rosettacode.org/wiki/Permutations#Fast_Lazy_Version  See also the good API of the Python combinations/permutations here: http://docs.python.org/3/library/itertools.html#itertools.permutations  note also the useful "r" and "repeat" arguments).

With such 7 functions/ranges you can do lot of things :-)

Bye,
bearophile

On Wednesday, 19 March 2014 at 01:29:16 UTC, bearophile wrote:
> - A function to compute the GCD on ulongs/longs/bigints is useful.
> (Issues 4125 and 7102).

Yeah, several methods work just fine if you change their declaration to isIntegral!T || is(typeof(T) == BigInt). gcd() is one of them.

Unfortunately, I don't trust rewriting isIntegral, so this sort of change would need to be on a function-by-function basis.

> - A function to compute the GCD on ulongs/longs/bigints is useful.
> (Issues 4125 and 7102).

A GCD is already in Phobos but it doesn't work with bigints. And putting it in the std.combinatorics is better.

Bye,
bearophile

Chris Williams:

> Yeah, several methods work just fine if you change their declaration to isIntegral!T || is(typeof(T) == BigInt). gcd() is one of them.
>
> Unfortunately, I don't trust rewriting isIntegral, so this sort of change would need to be on a function-by-function basis.

Don explained me that a GCD on BigInts should use a more efficient algorithm.

Bye,
bearophile

> nor to contain the most optimized algorithms around.

This D1 code adapted from C code is much more efficient (and in D2 with ranges and TypeTuple-foreach it could become more efficient and much shorter), but I think something like this is overkill for Phobos:

http://dpaste.dzfl.pl/cf97d15ade27

Bye,
bearophile