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 &lt; ceil(sqrt(MAX))){
        while (((++i)&lt;MAX) &amp;&amp; (!a[i]));

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

        for (unsigned long int k = i * i; k &lt; MAX; k += i)
            if (a[k])
                a[k] = false;
    }
    std::cout &lt;&lt; "done\n";
    for (unsigned long int i = 2; i &lt; MAX; i++)
        if (a[i])
            std::cout &lt;&lt; i &lt;&lt; "\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 &lt; ceil(sqrt(MAX))){

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

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

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

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

            if (a[k])
                //Open the following comment very carefully
                //std::cout &lt;&lt; "removes " &lt;&lt; 2 * k + 1 &lt;&lt; "\n";
                a[k] = false;
    }
    std::cout &lt;&lt; "done\n";
    for (unsigned long int i = 1; i &lt; MAX; i++)
        if (a[i])
            std::cout &lt;&lt; (2*i+1) &lt;&lt; "\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?
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 &lt; MAX; k += j)
    if (a[k])
        //Open the following comment very carefully
        //std::cout &lt;&lt; "removes " &lt;&lt; 2 * k + 1 &lt;&lt; "\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 &lt; MAX; k += j) {
    if (a[k]) {
        //no particular care is now required
        //std::cout &lt;&lt; "removes " &lt;&lt; 2 * k + 1 &lt;&lt; "\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 \$&lt;\sqrt{\text{MAX}}\$, but then there is this line:

while (((++i)&lt;MAX) &amp;&amp; (!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&lt;sqrtMAX &amp;&amp; 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 &lt; MAX; k += j) {


To this:

for (unsigned long int k = 2*i*(i+1); k &lt; 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 &lt;&lt; "done\n";
    for (i=1, j=2*i+1; i &lt; sqrtMAX; ++i, j+=2) {
        if (a[i]) {  // it's a prime
            std::cout &lt;&lt; j &lt;&lt; '\n';
            for (unsigned long int k = 2*i*(i+1); k &lt; MAX; k += j) {
                if (a[k]) {
                    a[k] = false;
                }
            }
        }
    }
    for ( ; i &lt; MAX; ++i, j+=2) {
        if (a[i]) { // it's a prime
            std::cout &lt;&lt; j &lt;&lt; '\n';
        }
    }
    free(a);
}


Further improvements

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


Use bits rather than bytes to store values.   
Reimplement with pointers rather than array indexing.   
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 &lt;typename T&gt;
class Bits 
{
public:
    Bits(size_t sz) :
        bitsperunit{8*sizeof(T)},
        units{sz/bitsperunit + 1},
        bits{new T[units]}
    {
        for (unsigned i=0; i &lt; units; ++i) {
            bits[i] = ~0;
        }
    }
    void clear(unsigned index) {
        unsigned n=index/bitsperunit;
        T mask=1u&lt;&lt;(index%bitsperunit);
        bits[n] &amp;= ~mask;
    }
    bool operator[](unsigned index) const {
        unsigned n=index/bitsperunit;
        T mask=1u&lt;&lt;(index%bitsperunit);
        return bits[n] &amp; 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&lt;uint8_t&gt; 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 &lt; sqrtMAX; ++i, j+=2) {
        if (a[i]) {  // it's a prime
            for (unsigned long int k = 2*i*(i+1); k &lt; MAX; k += j) {
                    a.clear(k);
            }
        }
    }
    std::cout &lt;&lt; "done\n";
    for (i=1, j=3 ; i &lt; MAX; ++i, j+=2) {
        if (a[i]) { // it's a prime
            std::cout &lt;&lt; j &lt;&lt; '\n';
        }
    }
}


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