Extending Sieve of Eratosthenes beyond a billion - follow-up

Follow up form Extending Sieve of Eratosthenes beyond a billion

Taking suggestions from the comments on the previous post:

Updated Code: This takes about 22 seconds for MAX as $10^9$ to reach "done". This is a significant improvement.

void sieve_of_eratosthenes(){
bool* a;
a = (bool*)malloc(MAX * sizeof(bool));
memset(a, true, MAX);
unsigned long int i = 1;
while (i < ceil(sqrt(MAX))){
while (((++i)<MAX) && (!a[i]));

if (2 * i >= MAX)//I am keeping this test. Can't Figure out if it is still relevant
break;

for (unsigned long int k = i * i; k < MAX; k += i)
if (a[k])
a[k] = false;
}
std::cout << "done\n";
for (unsigned long int i = 2; i < MAX; i++)
if (a[i])
std::cout << i << "\n";
free(a);
getchar();
getchar();
}

Further improved:

This takes about 24 seconds for MAX as $10^9$ to reach "done". This is a very strange behavior as I am processing only odd numbers now. So this should take less time(must be less than 22 seconds). Also, it prints from $3$ to $2*MAX + 1$ but even that isn't such a big margin. Good thing is with $10^9$ array size, I am able to find primes up to $2*10^9 + 1$.

void sieve_of_eratosthenes(){
bool* a;
a = (bool*)malloc(MAX * sizeof(bool));
memset(a, true, MAX);
unsigned long int i = 0;
while (i < ceil(sqrt(MAX))){

while (((++i)<MAX) && (!a[i]));

//if ((2 * i+1) >= MAX)//This test isn't going to be relevant anymore.
//break;
int j = 2 * i + 1;

//std::cout << "element is " << 2*i+1 << "\n";

for (unsigned long int k = i + j; k < MAX; k += j)

if (a[k])
//Open the following comment very carefully
//std::cout << "removes " << 2 * k + 1 << "\n";
a[k] = false;
}
std::cout << "done\n";
for (unsigned long int i = 1; i < MAX; i++)
if (a[i])
std::cout << (2*i+1) << "\n";

free(a);
getchar();
getchar();
}

Why is the strange behavior in time? Is there a limit to how much a malloc can allocate? Going beyond a billion for MAX still gives me error. How can I further improve this?

• Why didn't you follow the advice to use a std::vector<bool>? – πάντα ῥεῖ Dec 21 '15 at 10:03
• std::vector<bool> slowed me down very significantly. – piepi Dec 21 '15 at 10:04
• You may use reserve () or resize() once to prevent you from recurring reallocations. – πάντα ῥεῖ Dec 21 '15 at 10:06
• The standard specifies that std::vector<bool> has to use bit references. This mandatory inefficiency pretty much relegates its use to homework assignments and toy problems, since the overhead is difficult to eliminate even for current optimizing compilers. std::bitarray<> is better but fixed in size. Hence we have to roll our own or use third party libraries like boost. A simpler and more flexible alternative is templated bit accessor functions plus vector<unsigned>, as shown in Sieve of Eratosthenes - segmented with pre-sieving #2 – DarthGizka Dec 21 '15 at 12:39
• Oops, that would be std::bitset<>, not bitarray<> – DarthGizka Dec 21 '15 at 12:48

There are a few things you could do to improve your code.

Use braces

Constructions like this are a bug waiting to happen:

for (unsigned long int k = i + j; k < MAX; k += j)
if (a[k])
//Open the following comment very carefully
//std::cout << "removes " << 2 * k + 1 << "\n";
a[k] = false;

If anyone ever uncomments that line, the a[k] = false will then be outside the for loop and the program will not work correctly. Such errors are easily avoided by the use of braces:

for (unsigned long int k = i + j; k < MAX; k += j) {
if (a[k]) {
//no particular care is now required
//std::cout << "removes " << 2 * k + 1 << "\n";
a[k] = false;
}
}

Check for allocation failure

The code doesn't currently check for an allocation failure, but it should. Calls to malloc or calloc can fail, and your program should be robust enough not to fail if that happens.

Reduce the range of search for primes

We already know that the largest usable prime is $<\sqrt{\text{MAX}}$, but then there is this line:

while (((++i)<MAX) && (!a[i]));

We don't really need to search all the way to MAX. Also, for clarity, I'd rewrite that as a for loop:

// skip over composites to next prime
for (++i; i<sqrtMAX && a[i]; ++i)
{}

Skip more composite numbers

In the actual sieve, the code skips forward and then marks all odd multiples of the most recently discovered prime. However, the code currently doesn't skip as many as it could.

For instance, if the most recently discovered prime is 7, the code currently skips to 21, but we can do better. We can skip to 7*7 = 49 because any other smaller multiple of 7 will already have been eliminated. That transforms the inner loop from this:

for (unsigned long int k = i + j; k < MAX; k += j) {

To this:

for (unsigned long int k = 2*i*(i+1); k < MAX; k += j) {

Print as you go

Rather than searching the entire array at the end, only the part of the array after $\sqrt{\text{MAX}}$ really needs to be searched since all of the other primes were identified during sieveing. We can print them as they are discovered and save some time searching.

Increment rather than multiply

In the various loops, the expression 2*i+1 (stored in j in the first loop) is used to convert from an index to the actual numerical value it represents. It speeds things up a bit if the code simply increments j each loop increment.

Clean up the loop

Instead of using nested while loops, it seems to me to be cleaner if you use a for loop and inner if to determine whether the number is composite and should simply be skipped, or prime and should be sieved.

A new sieve:

Using all of these suggestions, I get this, which is about 8% faster on my machine than the original.

void sieve2()
{
bool *a = (bool*)malloc(MAX * sizeof(bool));
if (a == NULL) {
return;  // bad alloc
}
memset(a, true, MAX);
const unsigned long sqrtMAX = ceil(sqrt(MAX));
unsigned long i;
unsigned long j;
std::cout << "done\n";
for (i=1, j=2*i+1; i < sqrtMAX; ++i, j+=2) {
if (a[i]) {  // it's a prime
std::cout << j << '\n';
for (unsigned long int k = 2*i*(i+1); k < MAX; k += j) {
if (a[k]) {
a[k] = false;
}
}
}
}
for ( ; i < MAX; ++i, j+=2) {
if (a[i]) { // it's a prime
std::cout << j << '\n';
}
}
free(a);
}

Further improvements

There are many other things you might explore to further speed the code:

1. Use bits rather than bytes to store values.
2. Reimplement with pointers rather than array indexing.
3. Create a multithreaded version

Use classes

Since this is supposed to be C++, why not use a class? Here's a very simple class template that implements something like a bitset.

template <typename T>
class Bits
{
public:
Bits(size_t sz) :
bitsperunit{8*sizeof(T)},
units{sz/bitsperunit + 1},
bits{new T[units]}
{
for (unsigned i=0; i < units; ++i) {
bits[i] = ~0;
}
}
void clear(unsigned index) {
unsigned n=index/bitsperunit;
}
bool operator[](unsigned index) const {
unsigned n=index/bitsperunit;
return bits[n] & mask;
}
virtual ~Bits() {
delete[] bits;
}
private:
const unsigned bitsperunit;
const size_t units;
T *bits;
};

Now let's reimplement using that class:

void sieve3()
{
Bits<uint8_t> a(MAX);
const unsigned long sqrtMAX = std::ceil(std::sqrt(MAX));
unsigned long i;
unsigned long j;
for (i=1, j=2*i+1; i < sqrtMAX; ++i, j+=2) {
if (a[i]) {  // it's a prime
for (unsigned long int k = 2*i*(i+1); k < MAX; k += j) {
a.clear(k);
}
}
}
std::cout << "done\n";
for (i=1, j=3 ; i < MAX; ++i, j+=2) {
if (a[i]) { // it's a prime
std::cout << j << '\n';
}
}
}

On my machine this is still faster than either of the other two implementations (original or the one above).

• Perhaps worthy of a rethink: printing primes during the sieving effectively prevents testing with non-trivial values for MAX; given that the OP wants MAX to go beyond 10e9 that's a bit of a bummer. Who has the patience to wait for 50 million primes to be printed just in order to learn whether their code is actually correct? All the more so as printing via std::cout is slow as molasses. Love your avatar, though. ;-) – DarthGizka Dec 21 '15 at 18:16
• If I were writing this for my purposes, I wouldn't print them at all. However given a sieve of primes up to 2e9, measured speed on my computer showed printing-as-you-go is slightly faster. YMMV. – Edward Dec 21 '15 at 18:21

Your index math is all over the place, and you are calling std::sqrt() much more often than necessary, i.e. during every iteration of the outer loop instead of just once.

In order to clarify the issue it can be helpful to use a self-documenting type like uint32_t. Otherwise the code can become more complicated than necessary since the size of unsigned and unsigned long can vary from compiler to compiler. Also, the name MAX can clash with macros in commonly used headers; apart from that it is so broad that it doesn't really say what it means. For the sake of exposition I'll assume it to be a parameter called upper_limit.

uint32_t max_factor = uint32_t(std::sqrt(double(upper_limit)));

The cast to double is necessary for selecting the right overload. For upper limits beyond $2^{32}$ the precision of std::sqrt() (53 bits) can become insufficient and lead to wrong results, especially when sieving close to $2^{64}$. There's a separate topic for this very problem: Computing the square root of a 64-bit integer.

Also, you need to make up your mind whether you're doing a full sieve or an odds-only one. For odds-only sieves it is important to be very clear which values correspond to numbers and which ones correspond to bits; otherwise the code can become difficult to read and verify. One option is to use separate typedefs for numbers and sieve indices:

typedef uint32_t num_t;
typedef uint32_t idx_t;

The separation also allows it to change num_t to 64-bit while keeping indices as 32-bit, which would be the standard way for windowed/segmented sieves operating up to $2^{64}$.

In that way the computations can be expressed with greater clarity, for example like in this fragment (assuming an odds-only sieve where a set bit signifies 'composite' and which has been cleared already):

assert(upper_limit > 0);

idx_t max_sieve_bit = (upper_limit - 1) >> 1;

num_t max_factor = num_t(std::sqrt(double(upper_limit)));
idx_t min_factor_bit = 3 >> 1;     // an odds-only sieve doesn't have 2 in it
idx_t max_factor_bit = (max_factor - 1) >> 1;

for (idx_t bit = min_factor_bit; bit <= max_factor_bit; ++bit)
{
if (get_bit(sieve, bit))  // it's composite
continue;

num_t prime = (bit << 1) + 1;
idx_t start = (prime * prime) >> 1;
idx_t stride = (2 * prime) >> 1;

for (idx_t j = start; j <= max_sieve_bit; j += stride)
{
set_bit(sieve, j);
}
}

This version of the inner loop works only up to $2^{31}$ bits because otherwise j + stride can exceed UINT32_MAX even if max_sieve_bit does not. This means that the loop must be hardened against index math overflow if upper_limit can exceed UINT32_MAX.

get_bit() and set_bit() are just placeholders; you could implement them for your weapon of choice or make them macros or replace them with something else. It is strongly advisable to separate the sieve function from the output; otherwise it becomes impossible to write appropriate tests for the sieve function.

Also, don't fret over things like stride = (2 * prime) >> 1; modern compilers should optimise that to the moral equivalent of stride = prime automatically. The important thing is clarity, which should only be sacrificed when there is actual proof that it carries too heavy a price.

Index calculations and tests can be devilishly tricky in code like this, because it works in two different dimensions at the same time: numbers and bits that stand for odd numbers. Converting between the two offers lots of opportunities for fence-post errors, not least because the conversion from numbers to sieve bit indices is lossy. In my own code I tend to work almost exclusively in the realm of bits (sieve indexes) internally in order to simplify things, and only the external (exported, user-callable) functions take numbers as parameters.

Last but not least, setting a bit is lots faster than setting it only if it's not set (if (a[k]) a[k] = false; in your code). The reason is the conditional implied in the latter approach, and the fact that mis-predicted branches carry a heavy penalty in wasted CPU cycles.

More information on speeding up things can be found in Sieve of Eratosthenes - segmented to increase speed and range, which also has timings and links to compilable examples. Teaser: the linked article shows how to sieve the range up to $2^{32}$ in about 4.5 seconds; an improved version of the code - zrbj_sx_sieve32_v4.cpp - clocks 1.9 seconds with gcc and 2.3 seconds with VC++ 2013. The article also shows how to change the code for sieving up to $2^{64}$ and it explains why I laid out my sample code fragment the way I did: it makes it easier to add (and understand) the complications needed for high-speed sieving.

Short summary of techniques:

• packed odds-only bitmap (already doable with the fragment above)
• sieving in small segments that fit into the CPU's L1 cache
• remembering the working offset for each prime from segment to segment to avoid modulo division
• pre-sieving by small primes (blasting a fixed bit pattern over the sieve)

Techniques not shown which are employed by primesieve.org:

• using wheels for storage (extending the 'odds-only' idea to more small primes)
• bucketing for the working offsets

The last item is worthy of a more detailed explanation, because it is the reason why the primesieve.org program is the undisputed speed king.

Simple segmented sieves work well enough for small numbers but they start creaking loudly in the upper regions of the range. The reason is that sieving a segment requires them to go through all the prime factors up to the square root of the last number in the segment, regardless of how large or small the segment is. For the last segment just below $2^{64}$ that's more than 203 million primes, even if the segment holds just one number. The sieve does nothing with unneeded primes - it simply skips them after looking at the offset - but it does take quite a bit of time to skip hundreds of millions of them!

The trick of remembering offsets from segment to segment gets a lot more complicated. Instead of carrying millions of offsets from one segment to the next, each prime's working offset must be 'posted' to (stored for, associated with) the segment where the next multiple of that prime occurs.

The most important effect of this technique is that the sieving of any given segment can be restricted to primes that actually occur in that segment. On primesieve.org this technique is called 'bucketing' and it is the real battlefield for fast sieves beyond $2^{32}$. For example, the primesieve.org program can generate in excess of 1 GB of bucket/offset data internally.