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:%lun", 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=%dn", a);

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

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

            break;

    } 

    }

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

          printf("%d is primen", n);

    }

    gettimeofday(&edTime, NULL);

    printf("time: %u:%un", (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

Primality Testing — Part 3. A Further Optimized Method

The previous two posts cover the naive approaches for primality testing and a method with precomputation. This post introduces a new method which could perform better than previous approaches and it can be used together with the precomputation.

The optimization is based on the observation that all prime numbers are in the form of 6k+1 or 6k+5, with 2, 3 and 5 being the only exceptions. This is because all integers can be expressed in 6k+i for some integer k, and i  is from the set {0, 1, 2, 3, 4, 5}. And 2 divides 6k+0, 6k+2 and 6k+4, 3 divides 6k+3.

Therefore, the primality test for a given integer n can be done by checking if n is divisible by 2 and 3 first, and then checking if n is divisible for all numbers of form 6k+1 or 6k+5 and less than sqrt(n).

The idea of this optimization is actually similar to the Sieve of Eratosthenes algorithm described in previous post. We eliminated testing divisibility for numbers that is divisible by the prime factors.

In general, all prime numbers are in the form gk + i, where i < k and it represents the numbers that are coprime to g. In the case above, g = 6 and the coprime of 6 are 1 and 5. If i and g are not coprime, then i must be divisible by some prime factors of g.

We go through another example, suppose g = 2*3*5 = 30, then all integers can be expressed in the form 30*k + i where i is in the set {0, 1, 2….29}. Now we mark the numbers that divisible by 2, 3 and 5 respectively. 30*k + {0, 2, 4, 6, 8, 10, 12, …28} are divisible by 2; 30*k + {3, 6, 9, 12, 15, 18, 21, 24, 27} are divisible by 3; and 30*k + {5, 10, 15, 20, 25} are divisible by 5. The rest of the numbers {1, 7, 11, 13, 17, 19, 23, 29} are not divisible by 2, 3 and 5 and are therefore coprimes of 30.

So to test if a number n is prime number, we first test if the number is divisible by 2, 3, and 5 first, and then check if the number is divisible by numbers in the form 30*k + {1, 7, 11, 13, 17, 19, 23, 29}  up to sqrt(n).

Below is an implementation of the optimization for primality testing when g is set to 6,

/**

generate prime numbers

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int rdiv(int n, int m) {

    if (n%m ==0) {

        printf("%d is divisible by %dn", n, m);

        return 0;

    }

    return 1;

}

 

int prime(int n) {

    int i, k, ed, m;

    int cpr[2] = {1, 5};

    //check if we can divide the prime factors first

    if (n == 2 || n == 3) {return 1;}

    if (rdiv(n, 2) == 0) {

        return 0;

    } else if (rdiv(n, 3) == 0) {

        return 0;

    }

    //printf("ifprime: %d", ifprime);

    //check up to sqrt(n) of the form 6k + {coprime of 6} 

    ed = sqrt(n) + 1;

    for (k = 0; ;++k) {

        for (i = 0; i < 2; ++i) {

            m = 6*k + cpr[i];

            if (m < 3) {continue;}

            if (m > ed) {

                return 1;

            } 

            if (n%m == 0) {

                printf("%d is divisible by %dn", n, m);

                return 0;

            }

        }

    }

    return 1;

}

 

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

    int n = atoi(argv[1]);

    int i;

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

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

    } else {

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

    }

    gettimeofday(&edTime, NULL);

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

    return 0;

}

Save the code to prim4.c and compile the code with command below,

gcc -o prim4 prim4.c –lm

Below are some simple tests:

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

2 is prime number

time: 0:49

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 3

3 is prime number

time: 0:62

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 5

5 is prime number

time: 0:52

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2147483647

2147483647 is prime number

time: 0:267

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim4 2147483641

2147483641 is divisible by 2699

2147483641 is not prime number

time: 0:72

Note that the code can be further improved by adding the pre-computation we mentioned in part2.

Primality Testing — Part 2. Pre-computation with Sieve of Eratosthenes

Previous post covers two naive methods for primality testing. This post introduces an algorithm to generate prime numbers called Sieve of Eratosthenes. With this algorithm, an optimization can be applied to speed up the primality testing.

Sieve of Eratosthenes

Suppose we want to generate prime numbers up to N. The algorithm consists of two loops, and it can be described as below,

  • Create a mask for all integers up to N.
  • In the first loop, let x be 2, mark 2x, 3x, 4x for all n such that nx < N.
  • finding an integer bigger than x and not marked, let it be the new x, and repeat the step above
  • the first loop terminates when we cannot find an integer bigger than x and not marked
  • In the second loop, starting from 2, collects the numbers that are not marked, which are the prime numbers up to N.

Below is an implementation of the Sieve of Eratosthenes algorithm,

/**

generate prime numbers

**/

#include <stdio.h>

#include <stdlib.h>

 

#define MAXV 10000000

#define ININ 1000

 

unsigned char bits[MAXV/8];

 

void setbit(int pos) {

    bits[pos/8] |= (0x01 << (pos%8));

}

 

int testbit(int pos) {

    return (bits[pos/8] & (0x01 << (pos%8)));

}

 

void genprime(int n, int **pr, int *npr) {

   int i, j, k;

   for (i = 2; i < n; ++i) {

       if (!testbit(i)) {

           for (j = 2; i*j < n; ++j) {

               setbit(i*j);

           }

       }

   } 

   k = ININ;

   *pr = malloc(sizeof(int) * k);

   *npr = 0;

   for (i = 2, j = 0; i < n; ++i) {

       if (!testbit(i)) {

           (*pr)[j++] = i;

           (*npr)++;

           if (j == k) {

               k += ININ;

               *pr = realloc(*pr, sizeof(int)*k);

           }

       }

   }

}

 

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

    int n = atoi(argv[1]);

    int *pr;

    int npr;

    int i;

    struct timeval stTime, edTime;

    if (n > MAXV) {

        printf("the max n is %dn", MAXV);

        return 0;

    }

    gettimeofday(&stTime, NULL);

    genprime(n, &pr, &npr);

    printf("No. of prime numbers: %dn", npr);

    for (i = 0; i < npr; ++i) {

        printf("%d  ", pr[i]);

    }

    printf("n");

    gettimeofday(&edTime, NULL);

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

    if (pr != NULL) {

        free(pr);

    }

    return 0;

}

Note that the code above use a single bit to indicate if an integer is marked or not.

Save the code to prim2.c, and then compile the code with the command below,

gcc -o prim2 prim2.c

Below is a screenshot of some simple tests,


Figure 1. Simple Tests of Implementation for Sieve of Eratosthenes

Primality Testing with Sieve of Eratosthenes for Precomputation

One way to speed up primality testing is to precompute a list of prime numbers up to N. Then we first check if an input integer n is divisible by the prime numbers generated. If the number is not divisible by any of those prime numbers, we then proceed to the test mentioned in the naive approach post.

Below is an implementation of this precomputation approach,

/**

a primality test with precomputation

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

#define MAXV 10000000

#define ININ 1000

 

unsigned char bits[MAXV/8];

 

void setbit(int pos) {

    bits[pos/8] |= (0x01 << (pos%8));

}

 

int testbit(int pos) {

    return (bits[pos/8] & (0x01 << (pos%8)));

}

 

void genprime(int n, int **pr, int *npr) {

   int i, j, k;

   for (i = 2; i < n; ++i) {

       if (!testbit(i)) {

           for (j = 2; i*j < n; ++j) {

               setbit(i*j);

           }

       }

   } 

   k = ININ;

   *pr = malloc(sizeof(int) * k);

   *npr = 0;

   for (i = 2, j = 0; i < n; ++i) {

       if (!testbit(i)) {

           (*pr)[j++] = i;

           (*npr)++;

           if (j == k) {

               k += ININ;

               *pr = realloc(*pr, sizeof(int)*k);

           }

       }

   }

}

 

int prime(int n, int *pr, int npr, int x) {

    int i;

    int j;

    

    for (i = 0; i < npr; ++i) {

        if (n == pr[i]) {

            return 1;

        }

        if (n%pr[i] == 0) {

            printf("%d is divisble by %dn", n, pr[i]);

            return 0;

        }

    }

    if (n == 2) {

        return 1;

    }

    if (!(n&1)) {

        printf("%d is divisible by 2n", n);

        return 0;

    }

    j = sqrt(n);

    //to be safe, run to j+1

    i = (x+1) > 3 ? (x+1):3;

    if (!(i&1)) {

        ++i;   //make sure i starts as odd number

    }

    for (; i <= j + 1; i+=2) {

        if (n%i == 0) {

            printf("%d is divisible by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

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

    int m, n;

    int *pr;

    int npr;

    int i;

    struct timeval stTime, stTime2, edTime;

    if (argc < 3) {

        printf("./prim3 <primality test number> <prime number generation limit>n");

        return 0;

    }

    m = atoi(argv[1]);

    n = atoi(argv[2]);

    if (n > MAXV) {

        printf("the max n is %dn", MAXV);

        return 0;

    }

    gettimeofday(&stTime, NULL);

    genprime(n, &pr, &npr);

    gettimeofday(&stTime2, NULL);

    //printf("No. of prime numbers: %dn", npr);

/*    for (i = 0; i < npr; ++i) {

        printf("%dn", pr[i]);

    }*/

    if (prime(m, pr, npr, n)) {

        printf("%d is prime numbern", m);

    } else {

        printf("%d is not prime numbern", m);

    }

    gettimeofday(&edTime, NULL);

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

    printf("test prime time: %u:%un", (unsigned int)(edTime.tv_sec - stTime2.tv_sec), (unsigned int)(edTime.tv_usec - stTime2.tv_usec)); 

    if (pr != NULL) {

        free(pr);

    }

    return 0;

}

Save the source code to prim3.c, and then compile the code with the command below,

gcc -o prim3 prim3.c –lm

Below are some simple tests for the compiled program,

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

2 is prime number

total time: 0:144

test prime time: 0:48

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 3 10

3 is prime number

total time: 0:60

test prime time: 0:22

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483647 1000

2147483647 is prime number

total time: 0:469

test prime time: 0:302

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483643 1000

2147483643 is divisble by 3

2147483643 is not prime number

total time: 0:223

test prime time: 0:57

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim3 2147483641 1000

2147483641 is divisible by 2699

2147483641 is not prime number

total time: 0:250

test prime time: 0:80

References:
1. Wikipedia: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Primality Testing–Part 1. Naïve Approaches

Primality refers to determining if a given input number is prime or not. This post covers two naive approaches for primality testing.

Naive Approach 1

This most straightforward approach is to check if a given input integer n is divisible by integers from 2 to n-1. If a number is prime, then it should not divisible by any numbers in the range [2, n-1]. This is implemented as below,

/**

a naive method of determining if n is prime or not

divide from 2 to n-1

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int prime(int n) {

    int i;

    for (i = 2; i <= n-1; i++) {

        if (n%i == 0) {

            printf("%d is divisble by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

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

    int n = atoi(argv[1]);

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is primen", n);

    } else {

        printf("%d is not primen", n);

    }

    gettimeofday(&edTime, NULL);

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

}

Save the code to prim0.c. Compile the code with the command below,

gcc -o prim0 prim0.c –lm

Below are some sample testing results,

./prim0 2147483647

2147483647 is prime

time: 7:145569

 

./prim0 2147483640

2147483640 is divisble by 2

2147483640 is not prime

time: 0:46

 

./prim0 2

2 is prime

time: 0:51

 

./prim0 3

3 is prime

time: 0:49

Naive Approach 2: Improved Version of NA 1

The above algorithm can be improved. Firstly, if a number is divisible by an even integer, it must be divisible by 2. So we only need to check if the input number n is divisible by 2 and skip the rest of the even integers.

Secondly, if a number n is divisible by a number x, where x > sqrt(n). Then n = x*y, where y < sqrt(n), and y is also divisible by n. Therefore, we only need to check up to integer sqrt(n).

Below is the implementation with the two optimizations described,

/**

a naive method of determining if n is prime or not

divide from 2 to sqrt(n), skip the even numbers

**/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

 

int prime(int n) {

    int i, j;

    if (n == 2) {

        return 1;

    }

    if (!(n&1)) {

        //n is even

        printf("%d is divisible by 2n", n);

        return 0;

    }

    j = sqrt(n);

    //to be safe, run to j+1

    for (i = 3; i <= j + 1; i+=2) {

        if (n%i == 0) {

            printf("%d is divisible by %dn", n, i);

            return 0;

        }

    }

    return 1;

}

 

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

    int n = atoi(argv[1]);

    struct timeval stTime, edTime;

    gettimeofday(&stTime, NULL);

    if (prime(n)) {

        printf("%d is primen", n);

    } else {

        printf("%d is not primen", n);

    }

    gettimeofday(&edTime, NULL);

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

}

Note that sqrt is a floating point operation. We check until sqrt(n) + 1, just to be safe.

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

gcc -o prim1 prim1.c –lm

Below are some testing results,

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

2 is prime

time: 0:50

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 3

3 is prime

time: 0:60

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2147483647

2147483647 is prime

time: 0:308

 

roman10@ra-ubuntu-1:~/Desktop/integer/prim$ ./prim1 2147483641

2147483641 is divisible by 2699

2147483641 is not prime

time: 0:73

References:
1. Wikipedia Primality Test: http://en.wikipedia.org/wiki/Primality_test
2. Algorithms by Dasgupta, CH Papadimitriou and UV Vazirani