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).