Prime Number calculation bottleneck

Question

Inspired by another question here: Sum of Prime numbers under 2 million, I decided to try and implement Eratosthenes' Sieve from scratch using the algorithm described in the article.

unsigned int CalculateSumOfPrimes(const unsigned int number) {
    if(number <= 1) return 0;
    if(number == 2) return 2;
    std::vector<unsigned int> listOfPrimes;
    CalculateListOfPrimes(listOfPrimes, number);
    unsigned int sum = 0;
    for(unsigned int i = 0; i < listOfPrimes.size(); ++i) {
        sum += listOfPrimes.at(i);
    }
    return sum;
}

void CalculateListOfPrimes(std::vector<unsigned int>& container, const unsigned int number) {
    if(container.size() != 0) {
        std::cout << "Container must be empty!" << std::endl;
        return;
    }
    unsigned int current_prime_check = 2;
    std::vector<unsigned int> listOfNonPrimes;
    while(current_prime_check < number) {
        for(unsigned int i = current_prime_check; i < number; i += current_prime_check) {
            if(i == current_prime_check) {
                container.push_back(i);
                continue;
            }
            if(std::find(listOfNonPrimes.begin(), listOfNonPrimes.end(), i) != listOfNonPrimes.end()) continue;
            listOfNonPrimes.push_back(i);
        }
        ++current_prime_check;
        while(std::find(listOfNonPrimes.begin(), listOfNonPrimes.end(), current_prime_check) != listOfNonPrimes.end()) {
            ++current_prime_check;            
        }
    }

}

Currently the major bottleneck at large numbers (65536 or higher) is CalculateListOfPrimes(...) at 67% of the total 15 second run time.

Switching the usage from a std::vector to a std::list (since I'm only adding to the END of the the list, it seemed like a good idea), only saved 5%.

Looking at it, the inner while loop is going through EVERY possible number from the currently known non-prime number to the next unknown number, is there an algorithm to tell what the next number to check should be or whether or not the next number and potentially every number after it is guaranteed NOT to be in the list of primes?

Answer 1

Your main in-effeciency is here:

if(std::find(listOfNonPrimes.begin(), listOfNonPrimes.end(), i) != listOfNonPrimes.end()) 

For every non-prime you are a doing a search through the list of non primes. For both std::vector and std::list this is an O(n) operation.

The usual way of implementing this is to have an array where the removed elements are represented by their index into the array and thus represent a complexity of O(1) to both set and check if an element has already been set to false.

std::vector<bool> sieve(number + 1, true);

Then the following lines can be replaced:

if(std::find(listOfNonPrimes.begin(), listOfNonPrimes.end(), i) != listOfNonPrimes.end())
    continue;
listOfNonPrimes.push_back(i);

// Change too
sieve[i] = false;

Other things:

Every time through the loop you do a test:

if (i == current_prime_check) {
    container.push_back(i);
    continue;
}

This is only true the first element.

So hoist it out of the loop. Do the container.push_back(i); outside the loop and start i at the next position:

for (unsigned int i = 2 * current_prime_check; i < number; i += current_prime_check)
             //   ^^^^ Start one place up from where you were.

Replace the other find as well:

++current_prime_check;
while (std::find(listOfNonPrimes.begin(), listOfNonPrimes.end(), current_prime_check) != listOfNonPrimes.end()) {
    ++current_prime_check;
}

If you use the sieve this can be replaced by a simple test into the sieve:

++current_prime_check;
while (!sieve[current_prime_check]) {
    ++current_prime_check;
}

After Edit Comments:

We can make it quicker.

Your function that generates the prime:

void CalculateListOfPrimes(std::vector<unsigned int>& container, const unsigned int number)

You pass an array the function fills the array then you processes the array as the next stage. A better solution would be to pass a function that is applied to each prime as you find it. Then you can have a simple version that pushes the values into a container. Or alternatively you can have a function that just sums up the primes (so you do not even need the container or have to iterate over it).

So try:

template<typename Func>
void CalculateListOfPrimes(Func const& action, const unsigned int number);

Then you could use it like this:

struct PB
{
    std::vector<unsigned int>& cont;
    PB(std::vector<unsigned int>& c): cont(c) {}
    void operator()(unsigned int value) const { cont.push_back(value); }
};
…
std::vector<int> data;
CalculateListOfPrimes(PB(data), 2000000);

Or you can have a specific function

struct Sum
{
   unsigned int& sum;
   Sum(unsigned int& s) : sum(s) {}
   void operator()(unsigned int value) const { sum += value; }
};
…
unsigned int total;
CalculateListOfPrimes(Sum(total), 2000000);

Comments:

Don't use system("pause");. It's plainly platform specific. All you are doing is getting the program to wait. So read a line from the standard input.

std::getline(std::cin, line); // If you have other input you may need to flush first.

Hash define not a good move.

Prefer to use a function (there is no performance penalty as it will probably be inlined). If you must use a hash define then define it like this:

#define PAUSE  do { /* Action To Pause */ } while(false)

atoi()

primes = CalculateListOfPrimes(std::atoi(argv[2]));

Not standard. No way to check failures.

I would use boost::lexical_cast(argv[2]) it will throw an exception on failure.

Don't pass pointers around

DisplayResult(&primes, 0);

Change the function so you pass by const reference. This avoids the copy, stops you from accidentally modifying the array, and you can't accidentally pass a NULL so its always valid.

void DisplayResult(std::vector<bool> const& container, const unsigned int sum);

Truth values:

The result of the == operation is a boolean so you don't need to use a ternary operator on-top of that:

bool isSumArg = (std::strcmp("-s", argv[1]) == 0 ? true : false);

is exactly equivalent to the easier-to-read

bool isSumArg = (std::strcmp("-s", argv[1]) == 0);

Prime specific values:

A quick optimization here. You do not need this loop to go all the way to number.

You just need to loop as far as sqrt(number) after that nothing will affect the result.

while (current_prime_check < number) { … }

Make functions do one thing.

This function does lots of different things depending on the inputs that are not related.

void DisplayResult(std::vector<bool>* container, const unsigned int sum) {
    if(container == NULL && sum == 0) {
        return;
    }

    if(sum != 0) {
        std::cout << "Sum of Primes: " << sum << std::endl;
        return;
    } else if(container != NULL) {
        std::cout << "List of Primes:" << std::endl;
        for(unsigned int i = 1; i < container->size(); ++i) {
            if(container->at(i) == true) std::cout << std::setw(15) << std::right << i << std::endl;
        }
    }
}

Answer 2

You shouldn't append non-primes to a list and then perform a large number of find operations on it (which will all scan the entire container).

Instead a sieve is a table of whether a number is a multiple of some prime or not. Something like this which you then fill out.

std::vector<bool> sieve(limit + 1);

Then for example sieve[5] being false means 5 is a prime, whereas sieve[10], sieve[15] will be all set to true, as multiples of 5.