I implemented the Miller-Rabin primality test in Python/SAGE. My implementation is up at my GitHub.

Given an integer n, the Miller-Rabin primality test determines whether n is not prime or is likely prime. Here is how the test goes:

(1) Find integers k and q such that (2 ^ k)q = n – 1, where k > 0 and q is odd. Here is how you do that:

(1.1) Let k be a running count of the number of divisions that you make.

(1.2) Divide n – 1 by 2, and increment k. Let q be the quotient.

(1.2) While q is even, divide q by 2 and increment k.

For example, let’s say that n = 13. So, we divide 12 by 2. The quotient, of course, is 6. Because 6 is even, we divide 6 by 2, and we get 3. Because 3 is odd, we stop dividing, and q = 3. Because we divided two times, k = 2.

We can test whether these values of k and q satisfy the equation (2 ^ k)q = n – 1 by substituting in those values. And we do find that (2 ^ 2) * 3 = 12.

(2) Pick a random integer between 1 and n – 1. Call it r.

(3) Test whether (r ^ q) % n = 1.

(3.1) If the test is passed, then n is likely prime.

(3.2) If the test is not passed, then, for each integer j between 0 and k – 1, inclusive, test whether (r ^ ((2 ^ j) * q)) % n = n – 1.

(3.2.1) If any j passes that test, then n is likely prime.

(3.2.2) If no j passes that test, then n is not prime.

So, for example, again let n = 13, q = 3, and k = 2. Now let r = 9. Does (r ^ q) % n = 1?

Yes. It is true that (7 ^ 3) % 13 = 1. So, the Miller-Rabin primality test indicates that 13 is likely prime. This is a great result, because we know that 13 is prime.

So, why does this algorithm work? Well, it is a fact of number theory that if n > 2 is prime, then, for all 1 < a < n – 1,

(a ^ q) % n = 1

or

(a ^ ((2 ^ j) * q)) % n = n – 1, for some 0 <= j < k.

Because one or the other condition must hold true for all 1 < a < n – 1, it must hold true for a random 1 < a < n – 1. So, if we pick a random integer between 1 and n – 1 and it meets one of the above conditions, then we have evidence that n is prime.

However, we do not have proof that n is prime. For one thing, an above condition must hold for all integers between 1 and n – 1. But, when we take a random integer, we only examine one such integer. Moreover, an above condition must hold for all prime numbers greater than 2, but an above condition can hold for a composite number, too.

Nevertheless, although we don’t know that n is prime if the random integer satisfies one of the conditions, we do know that n is not prime if the random integer does not satisfy one of the conditions, because an above condition must hold for all integers on that interval.