Primality testing of a number is perhaps the most common problem concerning number theory that topcoders deal with. A prime number is a natural number which has exactly two distinct natural number divisors: 1 and itself. Some basic algorithms and details regarding primality testing and factorization can be found here.
The problem of detecting whether a given number is a prime number has been studied extensively but nonetheless, it turns out that all the deterministic algorithms for this problem are too slow to be used in real life situations and the better ones amongst them are tedious to code. But, there are some probabilistic methods which are very fast and very easy to code. Moreover, the probability of getting a wrong result with these algorithms is so low that it can be neglected in normal situations.
This article discusses some of the popular probabilistic methods such as Fermat’s test, Rabin-Miller test, Solovay-Strassen test.
All the algorithms which we are going to discuss will require you to efficiently compute (ab)%c ( where a,b,c are non-negative integers ). A straightforward algorithm to do the task can be to iteratively multiply the result with ‘a’ and take the remainder with ‘c’ at each step.
However, as you can clearly see, this algorithm takes O(b) time and is not very useful in practice. We can do it in O( log(b) ) by using what is called as exponentiation by squaring. The idea is very simple:
This idea can be implemented very easily as shown below:
iterations b x y
0 107 1 7
1 53 7 4
2 26 1 7
3 13 1 4
4 6 4 7
5 3 4 4
6 1 7 7
7 0 4 4
Now b becomes 0 and the return value of the function is 4. Hence (7107)%9 = 4.
The above code could only work for a,b,c in the range of type “int” or the intermediate results will run out of the range of “long long”. To write a function for numbers up to 10^18, we need to compute (a*b)%c when computing a*b directly can grow larger than what a long long can handle. We can use a similar idea to do that:
Here is some code which uses the idea described above ( you can notice that its the same code as exponentiation, just changing a couple of lines ):
We could replace x=(x*y)%c with x = mulmod(x,y,c) and y = (y*y)%c with y = mulmod(y,y,c) in the original function for calculating (ab)%c. This function requires that 2*c should be in the range of long long. For numbers larger than this, we could write our own BigInt class ( java has an inbuilt one ) with addition, multiplication and modulus operations and use them.
This method for exponentiation could be further improved by using Montgomery Multiplication. Montgomery Multiplication algorithm is a quick method to compute (a*b)%c, but since it requires some pre-processing, it doesn’t help much if you are just going to compute one modular multiplication. But while doing exponentiation, we need to do the pre-processing for ‘c’ just once, that makes it a better choice if you are expecting very high speed. You can read about it at the links mentioned in the reference section.
Similar technique can be used to compute (ab)%c in O(n3 * log(b)), where a is a square matrix of size n x n. All we need to do in this case is manipulate all the operations as matrix operations. Matrix exponentiation is a very handy tool for your algorithm library and you can see problems involving this every now and then.
Fermat’s Little Theorem
According to Fermat’s Little Theorem if p is a prime number and a is a positive integer less than p, then
Given below is a simple function implementing Fermat’s primality test:
Though Fermat is highly accurate in practice there are certain composite numbers p known as Carmichael numbers for which all values of a<p for which gcd(a,p)=1, (a(p-1))%p = 1. If we apply Fermat’s test on a Carmichael number the probability of choosing an a such that gcd(a,p) != 1 is very low ( based on the nature of Carmichael numbers ), and in that case, the Fermat’s test will return a wrong result with very high probability. Although Carmichael numbers are very rare ( there are about 250,000 of them less than 1016 ), but that by no way means that the result you get is always correct. Someone could easily challenge you if you were to use Fermat’s test :). Out of the Carmichael numbers less than 1016, about 95% of them are divisible by primes < 1000. This suggests that apart from applying Fermat’s test, you may also test the number for divisibility with small prime numbers and this will further reduce the probability of failing. However, there are other improved primality tests which don’t have this flaw as Fermat’s. We will discuss some of them now.
Key Ideas and Concepts
Fermat’s Little Theorem.
If p is prime and x2 = 1 ( mod p ), then x = +1 or -1 ( mod p ). We could prove this as follows:
Now if p does not divide both (x-1) and (x+1) and it divides their product, then it cannot be a prime, which is a contradiction. Hence, p will either divide (x-1) or it will divide (x+1), so x = +1 or -1 ( mod p ).
Let us assume that p – 1 = 2d * s where s is odd and d >= 0. If p is prime, then either as = 1 ( mod p ) as in this case, repeated squaring from as will always yield 1, so (a(p-1))%p will be 1; or a(s*(2r)) = -1 ( mod p ) for some r such that 0 <= r < d, as repeated squaring from it will always yield 1 and finally a(p-1) = 1 ( mod p ). If none of these hold true, a(p-1) will not be 1 for any prime number a ( otherwise there will be a contradiction with fact #2 ).
Let p be the given number which we have to test for primality. First we rewrite p-1 as (2d)s. Now we pick some a in range [1,n-1] and then check whether as = 1 ( mod p ) or a(s(2r)) = -1 ( mod p ). If both of them fail, then p is definitely composite. Otherwise p is probably prime. We can choose another a and repeat the same test. We can stop after some fixed number of iterations and claim that either p is definitely composite, or it is probably prime.
A small procedure realizing the above algorithm is given below:
Key Ideas and Concepts
Legendre Symbol: This symbol is defined for a pair of integers a and p such that p is prime. It is denoted by (a/p) and calculated as:
It is proved by Euler that:
(a/p) = (a((p-1)/2)) % p
So we can also say that:
(ab/p) = (ab((p-1)/2)) % p = (a((p-1)/2))%p * (b((p-1)/2))%p = (a/p)*(b/p)
Jacobian Symbol: This symbol is a generalization of Legendre Symbol as it does not require ‘p’ to be prime. Let a and n be two positive integers, and n = p1k1 * .. * pnkn, then Jacobian symbol is defined as:
So you can see that if n is prime, the Jacobian symbol and Legendre symbol are equal.
There are some properties of these symbols which we can exploit to quickly calculate them:
(a/n) = 0 if gcd(a,n) != 1, Hence (0/n) = 0. This is because if gcd(a,n) != 1, then there must be some prime pi such that pi divides both a and n. In that case (a/pi) = 0 [ by definition of Legendre Symbol ].
(ab/n) = (a/n) * (b/n). It can be easily derived from the fact (ab/p) = (a/p)(b/p) ( here (a/p) is the Legendry Symbol ).
if a is even, than (a/n) = (2/n)*((a/2)/n). It can be shown that:
(a/n) = (n/a)*(-1((a-1)(n-1)/4)) if a and n are both odd.
The algorithm for the test is really simple. We can pick up a random a and compute (a/n). If n is a prime then (a/n) should be equal to (a((n-1)/2))%n [ as proved by Euler ]. If they are not equal then n is composite, and we can stop. Otherwise we can choose more random values for a and repeat the test. We can declare n to be probably prime after some iterations.
Note that we are not interested in calculating Jacobi Symbol (a/n) if n is an even integer because we can trivially see that n isn’t prime, except 2 of course.
Let us write a little code to compute Jacobian Symbol (a/n) and for the test:
The various routines provided in the article can be highly optimized just by using bitwise operators instead of them. For example /= 2 can be replaced by “>>= 1″, “%2″ can be replaced by “&1″ and “*= 2″ can be replaced by “<<=1″. Inline Assembly can also be used to optimize them further.
Problems involving non-deterministic primality tests are not very suitable for the SRM format. But problems involving modular exponentiation and matrix exponentiation are common. Here are some of the problems where you can apply the methods studied above: