Miller-Rabin Primality Test
Problem: Determine if a number is prime, with an acceptably small error rate.
Solution: (in Python)
import random def decompose(n): exponentOfTwo = 0 while n % 2 == 0: n = n/2 exponentOfTwo += 1 return exponentOfTwo, n def isWitness(possibleWitness, p, exponent, remainder): possibleWitness = pow(possibleWitness, remainder, p) if possibleWitness == 1 or possibleWitness == p - 1: return False for _ in range(exponent): possibleWitness = pow(possibleWitness, 2, p) if possibleWitness == p - 1: return False return True def probablyPrime(p, accuracy=100): if p == 2 or p == 3: return True if p < 2: return False exponent, remainder = decompose(p - 1) for _ in range(accuracy): possibleWitness = random.randint(2, p - 2) if isWitness(possibleWitness, p, exponent, remainder): return False return True
Discussion: This algorithm is known as the Miller-Rabin primality test, and it was a very important breakthrough in the study of probabilistic algorithms.
Efficiently testing whether a number is prime is a crucial problem in cryptography, because the security of many cryptosystems depends on the use of large randomly chosen primes. Indeed, we’ve seen one on this blog already which is in widespread use: RSA. Randomized algorithms also have quite useful applications in general, because it’s often that a solution which is correct with probability, say, $ 2^$ is good enough for practice.
But from a theoretical and historical perspective, primality testing lied at the center of a huge problem in complexity theory. In particular, it is unknown whether algorithms which have access to randomness and can output probably correct answers are more powerful than those that don’t. The use of randomness in algorithms comes in a number of formalizations, the most prominent of which is called BPP (Bounded-error Probabilistic Polynomial time). The Miller-Rabin algorithm shows that primality testing is in BPP. On the other hand, algorithms solvable in polynomial time without randomness are in a class called P.
For a long time (from 1975 to 2002), it was unknown whether primality testing was in P or not. There are very few remaining important problems that have BPP algorithms but are not known to be in P. Polynomial identity testing is the main example, and until 2002 primality testing shared its title. Now primality has a known polynomial-time algorithm. One might argue that (in theory) the Miller-Rabin test is now useless, but it’s still a nice example of a nontrivial BPP algorithm.
The algorithm relies on the following theorem:
Theorem: if $ p$ is a prime, let $ s$ be the maximal power of 2 dividing $ p-1$, so that $ p-1 = 2^d$ and $ d$ is odd. Then for any $ 1 \leq n \leq p-1$, one of two things happens:
The algorithm then simply operates as follows: pick nonzero $ n$ at random until both of the above conditions fail. Such an $ n$ is called a witness for the fact that $ p$ is a composite. If $ p$ is not a prime, then there is at least a 3/4 chance that a randomly chosen $ n$ will be a witness.
We leave the proof of the theorem as an exercise. Start with the fact that $ a^ = 1 \mod p$ (this is Fermat’s Little Theorem). Then use induction to take square roots (the result has to be +/-1 mod p), and continue until you get to $ a^=1 \mod p$.
The Python code above uses Python’s built in modular exponentiation function pow to do fast modular exponents. The isWitness function first checks $ n^d = 1 \mod p$ and then all powers $ n^ = -1 \mod p$. The probablyPrime function then simply generates random potential witnesses and checks them via the previous function. The output of the function is True if and only if all of the needed modular equivalences hold for all witnesses inspected. The choice of endpoints being 2 and $ p-2$ are because 1 and $ p-1$ will always have exponents 1 mod $ p$.
Share this:
An introduction to Miller-Rabin algorithm and its Python
This article explains the Miller-Rabin primality test in cryptography. It consists of three parts. The first part gives the math background for this algorithm and adaptations to make it practical to real world use. The second part gives a python impeletion.
2. Math background for the Miller-Rabin primality test
The main reference for this section is An Introduction to Cryptography
1. Fermat’s primality test
First let us recall the Fermat’s little Theorem.
Let p be a prime number, then for any in unit of .
- For a general integer , is called a Fermat witness (for the compositeness) of if . is called a Fermat non-witness of if . Equivalently, is a called pseudoprime with Fermat non-witness .
- Let be the set of Fermat witness of , is called Carmichael number if is empty and n is composite.
It can be seen from Definition 2 clearly that if is a pseudo prime number with every unit being a Fermat non-witness, then is either a prime number or Carmichael number.
Fermat’s little Theorem leads to an algorithm called Fermat’s primality test.
we have the following naive deterministic algorithm.
def fun(n): for a in set1. n-1> if a^(n-1)!=1 mod n: return false return True The first problem of this algorithm is even if fun(n) is true, is either a prime number or a Carmichael number. As shown in here , there is infinite number of Carmichael numbers. Prime numbers and Carmichael number can not be distinguished simply.
The second problem is the time complexity.
The Fermat’s primality test has time complexity .
In fact, the first comes from the times of loops. To do the modular exponential . One may apply exponentiating by squaring algorithm to reduce the task to compute product of two number less than , this gives us a factor . The product of two numbers with digits less than can be treated by Montgomery modular multiplication with time complexity . So the time complexity of the loop body is .
The first factor gives a lot restriction to the time complexity and Fermat primality test an actual exponential time algorithm relative to .
One may attempt to modify the Fermat’s primality test algorithm to be a probabilistic one selecting different numbers in randomly for k times. This algorithm has time complexity .
Let . Let C be event that is a composite, be event that after the output running the “probabilistic” Fermat’s primality test is True , then Since can be Carmelson number, the upper bound of can be 1. Hence this algorithm will is not probabilistic.
2. Miller-Rabin primality test
Miller-Rabin test is an probabilistic polynomial algorithm. It improves the Fermat’s primality test at the two problems mentioned in previous section. It relied on the following proposition.
Let be a prime number. has exactly two solutions , for in .
Based on this we have the following proposition for prime numbers.
Let p be a prime number, . Suppose , where is an odd number, t is an integer. Then
One may prove it by considering the minimal element of form satisfying .
This is leads to following definition analoging Fermat’s witness.
For a general integer , is called a Miller-Rabin witness (for the compositeness) of . Let if
Otherwise is called a Miller-Rabin non-witness of .
Correspondingly the following algorithm is designed to test if a number is a Miller-Rabin non-witness of . We may assume .
def isMillerRabinNonWitness(a,n): find s,t for n-1 u=a^s mod n if u==1 or -1 mod n: return False loop t-1 times u=u*u mod n if u=-1 mod n: return False return True The time complexity of checking being Miller-Rabin non-witness is .
- Time complexity to get is .
- Time complexity to get is . [A similar analysis appears in Propersition 3]
- Time complexity for the loop is . [loop times, the loop body has time complexity ]
The Miller-Rabin primality probabilistic algorithm is designed as the following.
def fun(n): loop k times randomly choose different a in 1. n-1> if not isMillerRabinNonWitness(a,n): return False return TrueThis is a probabilistic polynomial algorithm of time complexity with respect .
we use the concept of precision in data science. The precision of Miller-Rabin primality probabilistic algorithm.
Let be event that is a composite, be event that is a prime. be event that after the output running Miller-Rabin primality probabilistic algorithm is True . The precision of Miller-Rabin primality probabilistic algorithm is defined as .
The Miller-Rabin primality probabilistic algorithm has good estimation for precision.
Let . L, then
Not like the case discussed in the remark 1.1, has the following description.
Let n be an integer, then , where is the Euler’s phi function.
Check here Theorem 9.11 for the proof.
Pick one estimation from estimations based on prime number theorem.
Let be the prime-counting function that gives the number of primes less than or equal to x, then
Check here for the proof.
By the Bayes formula, for a bit number ,let , we have
The last inequality require , the second last inequality used the theorem above.
Though the Miller-Rabin primality probabilistic algorithm can not be used for determine a number being prime theoretically. As noted here , the probability for a CPU has one error has a lower bound . So one may choose to guarantee the precision of the algorithm is suitable for practical use.
It was show here that admitting the generalized Riemann hypothesis, if is a Miller-Rabin witness of , then . This leads to a deterministic Miller-Rabin primality algorithm with polynomial time complexity ， which can be used to verify primality in theory.
3. Python implementation of Miller-Rabin primality test
In the implementation of the algorithm, I customized a unit test module to verify my coding. Generally it will verify if each function in the algorithm will run correctly by testing their result on some customized data.
We use the python code at the end of this section to implement the Miller Rabin test.
During the process of coding, it worth to note that the following code
u=pow(a,s,n) instead of u=a**s%n give the improvement of running speed.
For the algorithm using u=a**s%n it took 0.001s to get a 19 bits prime. For algorithm using u=pow(a,s,n) it took 4s to get a 19 bits prime.
This is because though a**s is implemented in python with times multiplication but the digit of increase to digits, which has time complexity with FFT. pow(a,s,n) is implemented with a combination of Montgomery modular multiplication and exponential by square and has a prefered time complexity .
It tooks 1.3 seconds for generate_prime() to get a prime of 1024 bits and 27 seconds for generate_prime() to get a prime of 2048 bits.
from random import randrange,getrandbits from math import log,ceil def nbit_odd(n): if n==1: return 1 return randrange(2**(n-1)+1,2**(n),2) def is_MRnonwitness(a,n): if a==0: return True s=n-1 t=0 while s&1==0: s=s>>1 t+=1 u=pow(a,s,n) if u-1==0 or u+1==n: return True for i in range(t-1): u=pow(u,2,n) if u+1==n: return True return False def generate_prime(n): # assert n>1 assert n>1 k=ceil(43*log(log(n))) while True: m=nbit_odd(n) find=True for i in range(k): a=getrandbits(n)%m if not is_MRnonwitness(a,m): find=False break if find: return m bnlucas / miller_rabin.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
| def miller_rabin ( n , k = 10 ): |
| if n == 2 : |
| return True |
| if not n & 1 : |
| return False |
| def check ( a , s , d , n ): |
| x = pow ( a , d , n ) |
| if x == 1 : |
| return True |
| for i in xrange ( s — 1 ): |
| if x == n — 1 : |
| return True |
| x = pow ( x , 2 , n ) |
| return x == n — 1 |
| s = 0 |
| d = n — 1 |
| while d % 2 == 0 : |
| d >>= 1 |
| s += 1 |
| for i in xrange ( k ): |
| a = randrange ( 2 , n — 1 ) |
| if not check ( a , s , d , n ): |
| return False |
| return True |
| # benchmark of 10000 iterations of miller_rabin(100**10-1); Which is not prime. |
| # 10000 calls, 11111 per second. |
| # 74800 function calls in 0.902 seconds |