General primality test in C++ 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 &lt;iostream&gt;
#include &lt;cstdlib&gt;
#include &lt;cmath&gt;

template &lt;class Integral&gt;
bool IsPrime (const Integral &amp;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 &amp;&amp; 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&lt;top &amp;&amp; is_prime; i+=6){
      is_prime = (n%(i+2) != 0) &amp;&amp; (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 &lt;&lt; "Format: " &lt;&lt; argv[0] &lt;&lt; " &lt;integer&gt;" &lt;&lt; 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).
Simple errors
This code includes &lt;cmath&gt; and &lt;cstdlib&gt;, 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 &lt;math.h&gt; and &lt;stdlib.h&gt; 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&lt;bool&gt; is quite cheap on my hardware):
#include &lt;atomic&gt;
#include &lt;cmath&gt;
#include &lt;cstdlib&gt;

template &lt;class Integral&gt;
bool is_prime(const Integral &amp;n){
    // Get rid of base cases
    if (n &lt; 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&lt;bool&gt; is_prime = true;
#pragma omp parallel for
    for (Integral i = 3;  i &lt; 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 &lt;iostream&gt;
int main(int argc, char* argv[]){
    // Check passed arguments
    if (argc != 2){
        std::cerr &lt;&lt; &quot;Format: &quot; &lt;&lt; argv[0] &lt;&lt; &quot; &lt;integer&gt;&quot; &lt;&lt; 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 &lt;atomic&gt;
#include &lt;cmath&gt;
#include &lt;cstdlib&gt;
#include &lt;initializer_list&gt;

template&lt;class Integral&gt;
bool is_prime(const Integral &amp;n)
{
    // Values under 30
    if (n &lt; 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&lt;bool&gt; is_prime = true;
#pragma omp parallel for
    for (Integral i = 30;  i &lt; 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).