A very quick (segmented) sieve of Eratosthenes.
Takes ~6s on a midrange machine to hit all 50 847 534 primes less than 1 billion, ending with 999999937
If you want to actually do anything with every prime (beyond counting them), there are three places to add a statement doing whatever is necessary with "lastP" -- one at the top (handles the special case '2'), one in the middle (handles "small" primes which are actually used to sieve), and one at the bottom (handles "large" primes which merely survive sieving).
In principle one can use a function object as parameter to allow generic operations on the primes, so add that if you want more general-purpose code (perhaps I'll do that later)
For higher limits you need to switch to wider types and follow the commented guidelines for the constants. For a fixed limit, changing B_SIZE may affect performance so if needed tune it (profile as you go, of course!). But this will get quite slow if you go to much higher numbers.
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 | #include <iostream>
#include <vector>
using namespace std;
typedef unsigned int UI;
const UI MAX_N = 1000000000;
const UI RT_MAX_N = 32000; // square of max prime under this limit should exceed MAX_N
const UI B_SIZE = 20000; // not sure what optimal value for this is;
// currently must avoid overflow when squared
// assumes p, b odd and p*p won't overflow!
void crossOut(UI p, UI b, bool mark[]) {
UI si = (p - (b % p)) % p;
if (si & 1)
si += p;
if (p*p > b)
si = max(si, p*p-b);
for (UI i = si/2; i < B_SIZE/2; i += p)
mark[i] = true;
}
// 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]
void calcPrimes() {
int pCount = 1; UI lastP = 2;
// *** do something with first prime (2) here ***
vector<UI> smallP; smallP.push_back(2);
bool mark[B_SIZE/2] = {false};
for (UI i = 1; i < B_SIZE/2; ++i)
if (!mark[i]) {
++pCount; lastP = 2*i+1;
// *** do something with the newest prime lastP here ***
smallP.push_back(lastP);
if (lastP*lastP <= B_SIZE)
crossOut(2*i + 1, 1, mark);
}
else mark[i] = false;
for (UI b = 1+B_SIZE; b < MAX_N; b += B_SIZE) {
for (UI i = 1; smallP[i]*smallP[i] < b+B_SIZE; ++i)
crossOut(smallP[i], b, mark);
for (UI i = 0; i < B_SIZE/2; ++i)
if (!mark[i]) {
++pCount; lastP = 2*i + b;
// *** do something with the newest prime lastP here ***
if (lastP <= RT_MAX_N)
smallP.push_back(lastP);
}
else mark[i] = false;
}
cout << "Found " << pCount << " primes in total.\n";
cout << "Recorded " << smallP.size() << " small primes, up to " << RT_MAX_N << '\n';
cout << "Last prime found was " << lastP << '\n';
}
int main() {
calcPrimes();
}
|
Made this a while ago, thought I'd share since it's faster than what's here. It is probably mainly of educational interest because for most applications, you only need to generate a list of primes once and can then use it forever (and such lists are easy to find for download if you don't want to generate them).
Should have mentioned: a big advantage of the segmented sieve is the memory usage -- if you don't need to keep the list of all primes around (i.e. can process them online, like for counting as illustrated above) then this method only keeps a boolean array of size
B_SIZEand a vector of small primes (those less thanRT_MAX_N), which for the range of practical numbers is very affordable.