Miller-Rabin Primality Test
Why I think this is cool:
- Very fast (much faster than naively checking all numbers up to , for large )
- Short, simple code
Contents:
- Quick-and-dirty Python Code
- The Test
- Guaranteed small-number tests
- Backstory - how I came across this
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 that you’d like to check is prime or not, first compute the numbers and such that (in other words, factor out all the powers of 2 from ). Then, for various different choices of , check that:
If, for any choice of , all of these congruences are false, then is definitely not prime.
If, for every tested, any the congruences are true, then is probably prime.
Furthermore, the more random values for we try, the more “probable” it is that is prime (but note that, unless we all the possible choices of according to the Miller test or the test cases described next, we can’t be certain is prime).
This pair of contrapositives may be easier to understand:
is prime at least one of the congruences is true for every
and
every equivalence is false for at least one is not prime
where .
Guaranteed small-number tests
Where this gets interesting is that people have pre-checked lots of numbers! Copy-pasting from Wikipedia:
When the number to be tested is small, trying all 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 < 2,047, it is enough to test = 2;
- if < 1,373,653, it is enough to test = 2 and 3;
- if < 9,080,191, it is enough to test = 31 and 73;
- if < 25,326,001, it is enough to test = 2, 3, and 5;
- if < 3,215,031,751, it is enough to test = 2, 3, 5, and 7;
- if < 4,759,123,141, it is enough to test = 2, 7, and 61;
- if < 1,122,004,669,633, it is enough to test = 2, 13, 23, and 1662803;
- if < 2,152,302,898,747, it is enough to test = 2, 3, 5, 7, and 11;
- if < 3,474,749,660,383, it is enough to test = 2, 3, 5, 7, 11, and 13;
- if < 341,550,071,728,321, it is enough to test = 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 , then is prime if and only if at least one congruence holds for each of and ! In other words, you only need to check the congruences for 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 (). 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 , which is preferable if 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 :).