>> 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 <[email protected]>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 >
