We have already seen that one of the easiest and most efficient ways to generate a list of prime numbers is via the Sieve of Eratosthenes. What if we just need to test if a specific number is prime, though?

The Naive Algorithm Optimized

Say we want to test whether the number N is prime or not. The first approach that comes to mind is probably to try to dive N by all numbers from 2 up to N-1. If the division is perfect in any of those cases (i.e., no remainder) then N is not prime. Else if it can’t be perfectly divisible by any of the numbers it’s prime.

However there are two improvements we can make to this algorithm right away. First of all there’s no need to try to divide N by all even numbers. If N is not evenly divisible by 2 then it won’t be by any other even number. So we just test for 2, and then for all odd numbers.

Second, we don’t need to go up to N-1. In fact going up to the square root of N would suffice. Why? Because if N is not prime then it can be de-composed into at least two factors, and one of those factors must be smaller than or equal to the square of N. If this was not true (i.e., both factors were larger than the square of N) the resulting number would be larger the N itself.

The algorithm with both of those tweaks already applied is below:

#include <stdio.h>
#include <math.h>

int isPrime(int n){
    int i;

    if (n==2)
        return 1;

    if (n%2==0)
        return 0;

    for (i=3;i<=sqrt(n);i+=2)
        if (n%i==0)
            return 0;

    return 1;
}

int main(){
    int i;
    int count = 0;

    for (i=2;i<10000000;i++){
        if(isPrime(i))
            count++;
    }

    printf("count=%dn",count);

return 0;
}

It found 664579 primes among the first 10,000,000 integers, and it took 15.329 seconds to do it.

Pre-Computation with Binary Search

Another approach would be to build a list of primes first using the Sieve of Eratosthenes, and then to perform binary search in the resulting array whenever we wanted to check if a number if prime. Here’s the algorithm:

#include <stdio.h>
#include <stdlib.h>

#define LIMIT 10000000 /*size of integers array*/
#define PRIMES 700000 /*size of primes array*/

int bSearch(int n,int array[], int start, int end){
    int middle = (start+end)/2;

    if (start>end)
        return 0;

    if (n==array[middle])
        return 1;

    if (n>array[middle])
        return bSearch(n,array,middle+1, end);
    else
        return bSearch(n,array,start,middle-1);
}

int main(){
    int i,j;
    int *primes,*numbers;
    int count = 0;

    primes = malloc(sizeof(int)*PRIMES);
    numbers = malloc(sizeof(int)*LIMIT);

    /*fill the array with natural numbers*/
    for (i=0;i<LIMIT;i++){
        numbers[i]=i+2;
    }

    /*sieve the non-primes*/
    for (i=0;i<LIMIT;i++){
        if (numbers[i]!=-1){
            for (j=2*numbers[i]-2;j<LIMIT;j+=numbers[i])
                numbers[j]=-1;
        }
    }

    /*transfer the primes to their own array*/
    j = 0;
    for (i=0;i<LIMIT&&j<PRIMES;i++)
        if (numbers[i]!=-1)
            primes[j++] = numbers[i];

    for (i=2;i<10000000;i++)
        if (bSearch(i,primes,0,j-1))
            count++;

    printf("count=%dn",count);

return 0;
}

This algorithm found the same 664579 among the first 10,000,000 integers, but it 2.584 seconds. A huge improvement, but with a drawback: we can only use it if we know forehand the maximum size of the integers will be testing, else we won’t find it on the pre-built array of primes.

Fermat’s Primality Test

Another approach to verify whether a number is prime or not is to use a probabilistic test. That is, you carry out some tests and determine with some degree of certainty whether or not it’s prime. The simplest and most famous of those tests is called Fermat’s Primality Test, and it’s based on Fermat’s little theorem (Pierre de Fermat was a french lawyer and took mathematics as a hobby!). The theorem goes like this:

If p is a prime number, then for any integer a, a^p − a will be evenly divisible by p.

Using our knowledge of modular arithmetic we could express the theorem as:

a^p ≡ a mod p

The same theorem can be expressed as:

a^(p-1) ≡ 1 mod p

Most people tend to use the second variation in their algorithms, but either will work fine.

The basic idea is this: if you choose a random integer and apply Fermat’s theorem with it to the number you are testing and the result is different from 1 (I am using the second variation here) it means that the number is composite for sure. If you get a 1, on the other hand, the number may or may not be prime. Repeat the test several times, and every time you get a 1 the probability of that number being prime increases. You would need to perform some statistics analysis, but I believe that with around 18 tests you would already get close to 99,9% accuracy.

The problem with Fermat’s test is that in order to make it reliable you need to compromise its speed. On top of that it’s not completely reliable. There’s a class of numbers, called Carmichael numbers, that will be reported as primes by Fermat’s test when in reality they are composite. These numbers are rare, but if you need 100% reliability then this would be a problem.

Anyway below is a simple implementation of Fermat’s test:

#include <stdio.h>
#include <stdlib.h>

int modulo (int a, int b, int c){
    int result = 1;

    while (b>0){
        result *= a;
        result %= c;
        b--;
    }

    return result;
}

int isPrime(int n){
    int i;
    int x;

    if (n==2)
        return 1;

    for (i=0;i<10;i++){
        x = rand()%n;
        if (x==0||x==1)
            x++;
        if (modulo(x,n-1,n)!=1)
            return 0;
    }

    return 1;
}

int main(){
    int i;
    int count = 0;

    for (i=2;i<10000;i++){
        if(isPrime(i))
            count++;
    }

    printf("count=%dn",count);

return 0;
} 

The modulo function above is quite slow, but I wrote it that way to make clear how the algorithm works. If you want to speed up the algorithm you can use exponentiation by squaring, as in the code below:

int modulo (int a, int b, int c){
    int result = 1;

    while (b){
        if (b&1){
          result *= a;
          result %= c;
        }
        b >>=1 ;
        a *= a;
        a %= c;
    }

    return result;
}

Finally, remember that you might need to tweak the number of tests you are performing on each number (you control this on the for loop inside the isPrime() function) to make sure you won’t end up with false positives. The algorithm with the optimized module function took 8.2 seconds to test all numbers below 10,000,000. This is not as good as the approach with binary search, but it’s twice as fast as our naive algorithm, and it also can be applied to numbers of any size. Notice that I used 10 tests per number, though, and this resulted in around 30 false positives, maybe with some Carmichael numbers in the middle (still this is a reliability rate of 99.99995% more or less).

There are two more probabilistic tests that are worth mentioning: the Miller-Rabin test and Solovay-Strassen one, but I’ll talk about those in the future, as they are slightly more complicated.