Skip to content

Gerry's World

A glimpse into my life

Miller-Rabin Primality Test

January 28, 2022

Why I think this is cool:

Contents:

Quick-and-dirty Python Code

Reproduced below is my github gist for a super short copy-paste solution (that’s admittedly difficult to read/understand).

For a more readable and comprehensive version, check out my other gist.

Finally, also check out rosettacode.org for implementations in tons of other languages and tons of other algorithms! (I just discovered the site and it’s amazing!)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

"""
Simple Miller-Rabin Primality Test
Gauranteed correct up to 318,665,857,834,031,151,167,461
(which is >2^64 so works for any uint64)
Usage:
print(is_probably_prime(17))                     # True
print(is_probably_prime(12345678910987654321))   # True
print(is_probably_prime(46))                     # False
"""
_known_primes = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)

def _is_composite(n, a, d, s):
    return ((x := pow(a, d, n)) != 1) and (x != n - 1) and (not any(
        (x := (x**2 % n)) == n - 1 for _ in range(s - 1)))

def is_probably_prime(n):
    if n in _known_primes:
        return True
    if any((n % p) == 0 for p in _known_primes) or n in (0, 1):
        return False
    d, s = n - 1, 0
    while d % 2 == 0:
        d, s = d // 2, s + 1  # compute (d, s) s.t. n = d * 2^s + 1
    return all(not _is_composite(n, a, d, s) for a in _known_primes)

The Test

See Wikipedia link, especially the section on small sets of bases. Also definitely check out the Numberphile video starring Matt Parker!!!

Given a number \(n\) that you’d like to check is prime or not, first compute the numbers \(d\) and \(s\) such that  \(n - 1 = 2^s d\)   (in other words, factor out all the powers of 2 from \(n-1\)). Then, for various different choices of \(a\), check that:

\[\begin{cases} a^d \equiv 1 \pmod{n}\\ a^{2^rd} \equiv -1 \pmod{n}, \;\; \forall 0 \le r < s \end{cases}\]

If, for any choice of \(a < n\), all of these congruences are false, then \(n\) is definitely not prime.

If, for every \(a\) tested, any the congruences are true, then \(n\) is probably prime.

Furthermore, the more random values for \(a\) we try, the more “probable” it is that \(n\) is prime (but note that, unless we all the possible choices of \(a\) according to the Miller test or the test cases described next, we can’t be certain \(n\) is prime).

This pair of contrapositives may be easier to understand:

\(n\) is prime \(\implies\) at least one of the congruences is true for every \(\mathbf{a}\)

and

every equivalence is false for at least one \(\mathbf{a}\) \(\implies\) \(n\) is not prime

where \(0\le a < n\).

Guaranteed small-number tests

Where this gets interesting is that people have pre-checked lots of numbers! Copy-pasting from Wikipedia:

When the number \(n\) to be tested is small, trying all \(a < 2(\ln n)^2\) is not necessary, as much smaller sets of potential witnesses are known to suffice. For example, Pomerance, Selfridge, Wagstaff and Jaeschke have verified that

  • if \(n\) < 2,047, it is enough to test \(a\) = 2;
  • if \(n\) < 1,373,653, it is enough to test \(a\) = 2 and 3;
  • if \(n\) < 9,080,191, it is enough to test \(a\) = 31 and 73;
  • if \(n\) < 25,326,001, it is enough to test \(a\) = 2, 3, and 5;
  • if \(n\) < 3,215,031,751, it is enough to test \(a\) = 2, 3, 5, and 7;
  • if \(n\) < 4,759,123,141, it is enough to test \(a\) = 2, 7, and 61;
  • if \(n\) < 1,122,004,669,633, it is enough to test \(a\) = 2, 13, 23, and 1662803;
  • if \(n\) < 2,152,302,898,747, it is enough to test \(a\) = 2, 3, 5, 7, and 11;
  • if \(n\) < 3,474,749,660,383, it is enough to test \(a\) = 2, 3, 5, 7, 11, and 13;
  • if \(n\) < 341,550,071,728,321, it is enough to test \(a\) = 2, 3, 5, 7, 11, 13, and 17.

(For a longer list, check out OEIS)

For example, people have shown through exhaustive testing that, as long as \(n<1373653\), then \(n\) is prime if and only if at least one congruence holds for each of \(a=2\) and \(a=3\)! In other words, you only need to check the congruences for \(a=2, 3\) and you’re good up to 1373653! (exclamation points used as excited punctuation, not factorials)

My interpretation, which I (baselessly) believe to be somewhat unique, is that this is almost like a super efficient compression algorithm to store a database of prime numbers. Some supercomputers crunched away finding a bunch of prime numbers, and rather than storing them as a giant list, you can instead just store the database in the form of a few key numbers (\(a\)). Then, on our local computer, we don’t have to search a giant database to check if a number is prime, we just use the key numbers, using the Miller-Rabin test as the “decompression algorithm”. Just 7 numbers contain the primality information for 341,550,071,728,321 numbers! Incredible!

Both the code samples above use these small-number tests instead of using randomly generated values of \(a\), which is preferable if \(n\) is larger than the largest guaranteed value.

Backstory - how I came across this

I was working on a Project Euler problem for fun the other day and, as often occurs in some problems, I needed a quick function for testing whether or not a number is prime.

I lost my previous Project Euler programs to my old laptop (where I had various convenience codes handy) so I made a quick Google search to find a python function for testing if a number is prime. To my disappointment, almost all the results were either a Sieve of Eratosthenes or of this naive form:

def is_prime(n):
    return all(n % k != 0 for k in range(2, math.ceil(math.sqrt(n))))

How could this naive implementation be so common? I felt like there had to be some clever math hacks for more efficient primality checks, especially for the reasonably large numbers I was checking. Finally, I expected that there should be standard libraries for this but didn’t come across them in my brief search. :/

I had come across a couple posts mentioning the Miller-Rabin test but it looked pretty complicated. Finally I gave in and decided to read what this mysterious Miller-Rabin test was and I was elated to discover actually it’s really simple!!!

And even better, I vaguely remembered seeing the Numberphile video about this a couple months ago, so I was really excited about this :).