Prime number calculator using Sieve of Eratosthenes

Question

I've been writing this prime number calculator, and I've gotten it pretty quick in comparison to most I've seen knocking about. This one uses the Sieve of Eratosthenes approach, and I've optimised the code as much as possible with my knowledge. A more experienced mind might be able to make it better.

Machine: 2.4GHz Quad-Core i7 w/ 8GB RAM @ 1600MHz

Compiler: clang++ main.cpp -O3

Benchmarks:

Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 100

Calculated 25 prime numbers up to 100 in 2 clocks (0.000002 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 1000

Calculated 168 prime numbers up to 1000 in 4 clocks (0.000004 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 10000

Calculated 1229 prime numbers up to 10000 in 18 clocks (0.000018 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 100000

Calculated 9592 prime numbers up to 100000 in 237 clocks (0.000237 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 1000000

Calculated 78498 prime numbers up to 1000000 in 3232 clocks (0.003232 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 10000000

Calculated 664579 prime numbers up to 10000000 in 51620 clocks (0.051620 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 100000000

Calculated 5761455 prime numbers up to 100000000 in 918373 clocks (0.918373 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 1000000000

Calculated 50847534 prime numbers up to 1000000000 in 10978897 clocks (10.978897 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$ ./a.out 4000000000

Calculated 189961812 prime numbers up to 4000000000 in 53709395 clocks (53.709396 seconds).
Caelans-MacBook-Pro:Primer3 Caelan$

Source:

#include <iostream> // cout
#include <cmath> // sqrt
#include <ctime> // clock/CLOCKS_PER_SEC
#include <cstdlib> // malloc/free

using namespace std;

int main(int argc, const char * argv[]) {
    if(argc == 1) {
        cout << "Please enter a number." << "\n";
        return 1;
    }
    long n = atol(argv[1]);
    long i;
    long j;
    long k;
    long c;
    long sr;
    bool * a = (bool*)malloc((size_t)n * sizeof(bool));

    for(i = 2; i < n; i++) {
        a[i] = true;
    }

    clock_t t = clock();

    sr = sqrt(n);
    for(i = 2; i <= sr; i++) {
        if(a[i]) {
            for(k = 0, j = 0; j <= n; j = (i * i) + (k * i), k++) {
                a[j] = false;
            }
        }
    }

    t = clock() - t;

    c = 0;
    for(i = 2; i < n; i++) {
        if(a[i]) {
            //cout << i << " ";
            c++;
        }
    }

    cout << fixed << "\nCalculated " << c << " prime numbers up to " << n << " in " << t << " clocks (" << ((float)t) / CLOCKS_PER_SEC << " seconds).\n";

    free(a);

    return 0;
}

Answer 1

This is written like C code; there are lots of things you're doing that would be much nicer using the tools C++ gives you.

First, though, I'm not a fan of a global using namespace std, so I removed it.

You're using lots of longs; it's best to replace them with int64_ts for consistency. Perhaps some size_ts would work well, too, but it doesn't seem like an improvement to the API. You should also initialize them in the loops or as late as possible.

Your bool * a should be replaced with a std::vector<bool>; not only is it far easier but it's likely faster and takes less space.

You don't need to return 0.

Your

int64_t c = 0;
for(int64_t i = 2; i < n; i++) {
    if(a[i]) {
        //std::cout << i << " ";
        c++;
    }
}

can be just

int64_t c = std::count(std::begin(a), std::end(a), true);

using the algorithm standard library.

C++ compilers can tell that sqrt(n) is constant, so you don't need to assign it to a variable.

It's standard in C++ to do ++i over i++, although it won't matter in this case.

You can initialize your vector with

std::vector<bool> a(n, true);
a[0] = a[1] = false;

to avoid the first loop.

Your loop

for(int64_t k = 0, j = 0; j <= n; j = (i * i) + (k * i), ++k) {
    a[j] = false;
}

can be simplified to

for(int64_t j = i * i; j < n; j += i) {
    a[j] = false;
}

This does loop over slightly different values (it avoids j=0), but those values were likely unintended anyway.

For your timings, you should use chrono. This is long and verbose, but it's a so much better library that you should just deal with the hassle.

You validate against argc but not whether the input is a valid number. Try using istringstream to do that:

int64_t n;
if (!(std::istringstream(argv[1]) >> n)) {
    std::cout << "Number invalid." << "\n";
    return 1;   
}

Note that this allows input like " 6-afsadf" or " 34 2 " which parse as 6 and 34 respectively; a full check is more complicated.

You should probably separate the prime-counting code from the input-output code, although this requires a bit of hassle to return multiple values:

namespace ch = std::chrono;

std::pair<int64_t, ch::steady_clock::duration> primes_count_and_time(int64_t n) {
    ... // code

    auto elapsed = ch::steady_clock::now() - start;
    int64_t count = std::count(std::begin(a), std::end(a), true);
    return {count, elapsed};
}

int main(int argc, const char * argv[]) {
    ... // code

    int64_t count;
    ch::steady_clock::duration elapsed;
    std::tie(count, elapsed) = primes_count_and_time(n);

    ... // code
}

Now we have some simple optimizations to apply; you can start by dropping the odd values from the vector:

for(int64_t i = 1; i <= sqrt(n / 2 + 1); ++i) {
    if(a[i]) {
        // The start is:
        //      compress(uncompress(i)²)
        //          = ((2*i+1)² - 1) / 2
        //          = 2 * i * (i+1)
        for(int64_t j = 2*i*(i+1); j < n / 2 + 1; j += 2*i+1) {
            a[j] = false;
        }
    }
}

This takes 3.3 seconds for \$n = 10^9\$ for me, where the old version took 11.2. All of the changes are below:

#include <algorithm> // count
#include <chrono>    // steady_clock
#include <cmath>     // sqrt
#include <iostream>  // cout
#include <iterator>  // begin/end
#include <sstream>   // istringstream
#include <tuple>     // tee
#include <utility>   // pair
#include <vector>    // vector

namespace ch = std::chrono;

std::pair<int64_t, ch::steady_clock::duration> primes_count_and_time(int64_t n) {
    std::vector<bool> a(n / 2 + 1, true); // 2, 3, 5, 7... might be prime

    auto start = ch::steady_clock::now();

    for(int64_t i = 1; i <= std::sqrt(n / 2 + 1); ++i) {
        if(a[i]) {
            // The start is:
            //      compress(uncompress(i)²)
            //          = ((2*i+1)² - 1) / 2
            //          = 2 * i * (i+1)
            for(int64_t j = 2*i*(i+1); j < n / 2 + 1; j += 2*i+1) {
                a[j] = false;
            }
        }
    }

    auto elapsed = ch::steady_clock::now() - start;
    int64_t count = std::count(std::begin(a), std::end(a), true);
    return {count, elapsed};
}

int main(int argc, const char *argv[]) {
    if(argc == 1) {
        std::cout << "Please enter a number." << "\n";
        return 1;
    }

    int64_t n;
    if (!(std::istringstream(argv[1]) >> n)) {
        std::cout << "Number invalid." << "\n";
        return 1;   
    }

    int64_t count;
    ch::steady_clock::duration elapsed;
    std::tie(count, elapsed) = primes_count_and_time(n);

    double seconds = ch::duration_cast<ch::duration<double>>(elapsed).count();
    std::cout << std::fixed
              << "\nCalculated " << count
              << " prime numbers up to " << n
              << " in " << elapsed.count()
              << " clocks (" << seconds << " seconds).\n";
}