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

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

• 1. I didn't know that c libraries required std:: namespace, which I'm really glad to know. 2. The point of using templates was so any kind of Integral can be checked, but I forgot to do the same inside the function. 3. The return of 0 or 1 is for bash condition statements. 4. I've never heard about the -march= option and when to use it. Also, I usually don't use -O3 (always -O2). 5. I considered parallelizing the algorithm, but I wanted to go as far as I could using a single thread. 6. I thought of a different way of unrolling the loop, but yours works as well. – Francisco Gallego Salido Aug 4 '17 at 17:01
• I didn't mean to suggest that you must use -O3 - I noted that only to be clear how I tested. -march=native can be very useful - it tells the compiler to generate the best code it can for the processor on which we're compiling, but which might not run elsewhere. -mtune=native is similar, but the code will at least run on other versions of the processor (possibly not as quickly). – Toby Speight Aug 7 '17 at 9:02
• I wasn't sure whether you'd considered parallelising the the test, so I looked for improvements in single-threaded and parallel execution - I hope it was helpful to you. – Toby Speight Aug 7 '17 at 9:22
• I now assume that the -march option stands for "machine architecture", which sound really useful for personal use, I really appreciate that! Thank you very much for your answer, I didn't thought I would have learnt that much in a single post. – Francisco Gallego Salido Aug 7 '17 at 18:36

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.