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

## One comment on “Primality Testing — Part 4. Fermat’s Little Theorem and A Probabilistic Approach”

1. sir i want a program to check the primality condition if the given number is even

You may use these HTML tags and attributes: `<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> `