General primality test in C++

Question

I have this code, which I think can be optimized even more, but I can't think of a way of optimizing the main loop to make it faster and more comprehesive. Here's the code:

#include <iostream>
#include <cstdlib>
#include <cmath>

template <class Integral>
bool IsPrime (const Integral &n){
   // Get rid of base cases
   if (n==2 || n==3) return true;

   // Initialize the boolean to an intuitive correct value
   bool is_prime = n%2 != 0 && n%3 != 0;
   unsigned long top = (unsigned long) sqrt(n)+1;

   // This is the loop I want to optimize
   for (unsigned long i=3; i<top && is_prime; i+=6){
      is_prime = (n%(i+2) != 0) && (n%(i+4) != 0);  // kind of loop unrolling...
   }

   return is_prime;
}

// Main program
int main (int argc, char* argv[]){
   // Check passed arguments
   if (argc != 2){
      std::cerr << "Format: " << argv[0] << " <integer>" << std::endl;
      return -1;
   }

   if (IsPrime(strtoul(argv[1], NULL, 0))) return 0;
   else return 1;
}

As you can see, I could reduce the amount of iterations by six times, but going further was worthless.

For the biggest 64-bit unsigned prime number, 18446744073709551557, it spends about 15-17 seconds in an i5, 1.70 GHz processor (which I think can be optimized more).

Answer 1

Simple errors

This code includes <cmath> and <cstdlib>, but then calls unqualified sqrt and strtoul. Implementations are not required to import these names to the global namespace, so they should be qualified with std:: (or include <math.h> and <stdlib.h> instead).

An unsigned long may not be big enough to hold the integer part of sqrt(n) (for example, if Integer is a 128-bit type, and unsigned long is only 32 bits). The only safe way (short of heavy template metaprogramming) is to use Integer for the counter.

It's probably worth explicitly returning false for 0 and 1, and for any negative values that crazy people feel compelled to test.

Style

We usually use snake_case (or in some projects, camelCase) for functions, and reserve PascalCase for composite types. But go along with the conventions of your collaborators where necessary.

Performance

I did manage to get a small improvement (from 9½ seconds down to about 9 here, compiled with gcc -march=native -03) by splitting the condition into two (I think this might be allowing the compiler to execute both tests together).

The best results were obtained by simply parallelizing the algorithm (std::atomic<bool> is quite cheap on my hardware):

#include <atomic>
#include <cmath>
#include <cstdlib>

template <class Integral>
bool is_prime(const Integral &n){
    // Get rid of base cases
    if (n < 2) return false;
    if (n==2 || n==3) return true;

    if (n%2 == 0 || n%3 == 0) return false;

    auto const top = (Integral)std::sqrt(n) + 1;
    std::atomic<bool> is_prime = true;
#pragma omp parallel for
    for (Integral i = 3;  i < top;  i += 6) {
        if (!is_prime) continue;
        if (n%(i+2) == 0) is_prime = false;
        if (n%(i+4) == 0) is_prime = false;
    }
    return is_prime;
}

#include <iostream>
int main(int argc, char* argv[]){
    // Check passed arguments
    if (argc != 2){
        std::cerr << "Format: " << argv[0] << " <integer>" << std::endl;
        return -1;
    }

    if (is_prime(std::strtoul(argv[1], NULL, 0))) return 0;
    else return 1;
}

This got my time down to under 1½ seconds elapsed time, making full use of 8 cores.

I don't know how I can break instead of continue in an OpenMP parallel loop, but I found it didn't hurt too much with a simple composite number such as 625, which evaluated in around 20ms.

Loop unrolling

I did gain a speedup of a further 20% or so by unrolling the loop to consider 30 factors at a time:

#include <atomic>
#include <cmath>
#include <cstdlib>
#include <initializer_list>

template<class Integral>
bool is_prime(const Integral &n)
{
    // Values under 30
    if (n < 2) return false;
    for (Integral i: { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 })
        if (n%i == 0) return n==i;

    auto const top = (Integral)std::sqrt(n) + 1;
    std::atomic<bool> is_prime = true;
#pragma omp parallel for
    for (Integral i = 30;  i < top;  i += 30) {
        if (!is_prime) continue;
        if (n%(i+1) == 0) is_prime = false;
        if (n%(i+7) == 0) is_prime = false;
        if (n%(i+11) == 0) is_prime = false;
        if (n%(i+13) == 0) is_prime = false;
        if (n%(i+17) == 0) is_prime = false;
        if (n%(i+19) == 0) is_prime = false;
        if (n%(i+23) == 0) is_prime = false;
        if (n%(i+29) == 0) is_prime = false;
    }
    return is_prime;
}

int main()
{
    return is_prime(0)
        || is_prime(1)
        || !is_prime(2)
        || is_prime(49)
        || is_prime(961)
        || !is_prime(18446744073709551557u);
}

(I've also simplified main() here, to just run auto tests rather than requiring an argument).

Answer 2

To my eye, you can make some improvements when looking at base cases:

bool IsPrime (const Integral &n){

  // Get rid of base cases.
  if (n <= 1) {
    return false;
  }

  // Handle even numbers.
  if (n % 2 == 0) {
    return n == 2;
  }  

  // Handle multiples of 3.
  if (n % 3 == 0) {
    return n == 3;
  }

  // ...

That diverts all multiples of 2 and 3 before you have to set up your loop only not to run it: && is_prime. I also find it easier to read, though YMMV.

Your loop unrolling is a bit strange to my eye. You start at 3, which you have already checked for, instead of 5 and have to compensate in every repetition of your loop. Checking that && is_prime every time seems unnecessary as you can return false; as soon as any divisor is found.

for (unsigned long i = 5; i < top; i += 6) {
  if ((n % i != 0) || (n % (i + 2) != 0) {
    return false;
  }
}

// If we get here then the number is prime.
return true;

I haven't run any timings, but this has a simplified end-of-loop check and one less addition per loop so it might shave a little off the time.