Primality Testing — Part 4. Fermat’s Little Theorem and A Probabilistic Approach

Previous posts on primality testing are all trying to factor the input number. For a very large number with few hundreds of bits, it is difficult to compute in that manner.
This post introduce a probabilistic approach based on Fermat’s primality test.

Fermat’s Little Theorem

Fermat’s primality test is based on Fermat’s Little Theorem,

If p is prime, then for every 1 <= a < p,

[ a ^ ( p - 1) ] mod p = 1

Fermat’s Little Theorem suggests a method to test if a number is prime or not. For a big input number p, it is not desirable to test every a. Normally, we randomly pick up an a, and if it satisfy the equation above, we claim it is prime, otherwise, it is not.

However, Fermat’s little theorem is a necessary condition for a number to be prime, but not sufficient. In other words, every prime number satisfy fermat’s little theorem, but it is also possible to find an a for a composite number n where [ a ^ (p - 1) ] mod n = 1.

It can be proven that most values of a will fail the test for a composite number. And if we pick several random a and they all pass the test, the possibility of making a wrong decision is rare. This algorithm cannot guarantee that positive primality test is 100% correct, but it is fast the possibility of wrong is low enough for many practical use cases.

Implementation

Below is a simple C implementation based on the idea above,

#include <stdio.h>

#include <stdlib.h>

unsigned long modexp(int x, int y, int n) {

    unsigned long i;

    if (y == 0) {

        return 1%n;

    }

    if (y == 1) {

        return x%n;

    } 

    i = modexp(x, y/2, n);

    if (y&0x01) {    //odd

        //return (x*i*i)%n;

//        printf("%d:%lu:%lu\n", x, i, (x%n*i*i));

        return x%n*i%n*i%n;

    } else {         //even

        return i*i%n;

    }

}

int main(int argc, char **argv) {

    int n;

    int a;

    int i;

    n = atoi(argv[1]);

    srand(time(NULL));

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    for (i = 0; i < 20 && i < n; ++i) {

        a = rand()%100;          

        if (a == 0) {

            a = 1;

        }

        //printf("a=%d\n", a);

    if (modexp(a, n-1, n) != 1) {

            printf("%d is not prime, %d\n", n, a);

            break;

    } 

    }

    if (i == 20 || i == n) {

          printf("%d is prime\n", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%u\n", (unsigned int)(edTime.tv_sec - stTime.tv_sec), (unsigned int)(edTime.tv_usec - stTime.tv_usec)); 

}

Note that the code uses the modular exponentiation algorithm. Due to the integer overflow, the largest input the program can accept is 65536 = 2^16.

Save the code to prim5.c and compile it using the command below,

gcc -o prim5 prim5.c

Below are some simpe tests,

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 65521

65521 is prime

time: 0:69

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 2

2 is prime

time: 0:53

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 65536

65536 is not prime, 8

time: 0:53

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 63703

63703 is prime

time: 0:69

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 61493

61493 is prime

time: 0:67

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim5 13

13 is prime

time: 0:54

References:
1. Wikipedia Primality test: http://en.wikipedia.org/wiki/Primality_test
2. Algorithms