A simple sieve in Phobos?

Ziad Hatahet hatahet at gmail.com
Sun Mar 30 20:53:37 PDT 2014


>>      if (p ^^ 2 > b)
>>            si = si.max(p ^^ 2 - b);

Can't that be written as:

auto t = p ^^ 2 - b;
if (t > 0 && t > si) si = t;


--
Ziad


On Tue, Mar 18, 2014 at 7:23 AM, bearophile <bearophileHUGS at lycos.com>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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.puremagic.com/pipermail/digitalmars-d/attachments/20140330/5688cd34/attachment.html>


More information about the Digitalmars-d mailing list