An introduction to Miller-Rabin algorithm and its Python

1. Introduction

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.

Theorem 1.

(Fermat’s little Theorem).

Let p be a prime number, then $a^{p-1}=1$ for any $a\in (\bZ/ p\bZ)^*$ in unit of $\bZ/ p\bZ$ .

Definition 2.

$\textrm{ }$

For a general integer $n$ , $a\in (\bZ/ n\bZ)^*$ is called a Fermat witness (for the compositeness) of $n$ if $a^{n-1}\neq{1}$ . $a$ is called a Fermat non-witness of $n$ if $a^{n-1}=1$ . Equivalently, $n$ is a called pseudoprime with Fermat non-witness $a$ .
Let $F$ be the set of Fermat witness of $n$ , $n$ is called Carmichael number if $F$ is empty and n is composite.

It can be seen from Definition 2 clearly that if $n$ is a pseudo prime number with every unit being a Fermat non-witness, then $n$ 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 set{1,...,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, $n$ 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.

Propersition 3.

The Fermat’s primality test has time complexity $\tilde\sO(n)$ .

In fact, the first $n$ comes from the times of loops. To do the modular exponential $a^{n-1} \mod n$ . One may apply exponentiating by squaring algorithm to reduce the task to compute product of two number less than $n$ , this gives us a factor $\log n$ . The product of two numbers with digits less than $\log n$ can be treated by Montgomery modular multiplication with time complexity $\tilde\sO(\log n)$ . So the time complexity of the loop body is $\tilde\sO(\log^2 n)$ .

The first factor $n$ gives a lot restriction to the time complexity and Fermat primality test an actual exponential time algorithm relative to $\log n$ .

Remark 4.

One may attempt to modify the Fermat’s primality test algorithm to be a probabilistic one selecting different numbers in $\{1,..,,n-1\}$ randomly for k times. This algorithm has time complexity $\tilde\sO(k\cdot\log^2 n)$ .

Let $N_F=\{\textrm{Fermat non-Witness of } n\}$ . Let C be event that $n$ is a composite, $T_k$ be event that after the output running the “probabilistic” Fermat’s primality test is True, then $P(T_k|C)\leq(\frac{|N_F|}n)^k.$ Since $n$ can be Carmelson number, the upper bound of $P(T_k|C)$ can be 1. Hence this algorithm will is not probabilistic.<p></p>

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.

Propersition 5.

Let $p$ be a prime number. $x^2-1=0$ has exactly two solutions $1$ , $-1$ for $x$ in $\bZ/p\bZ$ .

Based on this we have the following proposition for prime numbers.

Propersition 6.

Let p be a prime number, $a\in (\bZ/ n \bZ)^*$ . Suppose $p-1=s\cdot2^t$ , where $s$ is an odd number, t is an integer. Then

$\begin{equation} \begin{split} &a^s=\pm1\textrm{ or} \\ &a^{s2^k}=-1 \textrm{ for some }k\in\{1,\dots,t-1 \} \end{split} \end{equation}$

One may prove it by considering the minimal element of form $a^{s2^k}$ satisfying $a^{s2^k}=1$ .

This is leads to following definition analoging Fermat’s witness.

Definition 7.

For a general integer $n$ , $a\in (\bZ/ n\bZ)^*$ is called a Miller-Rabin witness (for the compositeness) of $n$ . Let $n-1=s2^{t}$ if
$\begin{equation} \begin{split} a^s&=\pm1\textrm{ or} \\ a^{s2^k}&=-1 \textrm{ for some }k\in\{1,\dots,t-1 \} \end{split} \end{equation}$
Otherwise $a$ is called a Miller-Rabin non-witness of $n$ .

Correspondingly the following algorithm is designed to test if a number $a$ is a Miller-Rabin non-witness of $n$ . We may assume $a\in \{1,\dots ,n-1\}$ .

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 $a$ being Miller-Rabin non-witness is $\tilde\sO(\log ^2 n)$ .

Time complexity to get $s$ is $\sO(\log n)$ .
Time complexity to get $u$ is $\tilde\sO(\log^2 n)$ . [A similar analysis appears in Propersition 3]
Time complexity for the loop is $\tilde\sO(\log^2 n)$ . [loop $\sO(\log n)$ times, the loop body has time complexity $\tilde\sO(\log n)$ ]

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 True

This is a probabilistic polynomial algorithm of time complexity $\tilde\sO(k\cdot\log^2 n)$ with respect $\log n$ .

we use the concept of precision in data science. The precision of Miller-Rabin primality probabilistic algorithm.

Definition 8.

Let $C$ be event that $n$ is a composite, $P$ be event that $n$ is a prime. $T_k$ 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 $P(P|T_k)$ .

The Miller-Rabin primality probabilistic algorithm has good estimation for precision.
Let $N=\{\textrm{Miller-Rabin non-Witness of } n\}$ . L, then
$P(T_k|C)\leq(\frac{|N|}n)^k,$
and
$P(P|T_k)=1-P(T_k|C)\geq1-(\frac{|N|}n)^k.$

Not like the case discussed in the remark 1.1, $|N|$ has the following description.

Theorem 9.

Let n be an integer, then $|N|\leq \phi(n)/4<n/4$ , where $\phi(n)$ is the Euler’s phi function.

Check here Theorem 9.11 for the proof.

Pick one estimation from estimations based on prime number theorem.

Theorem 10.

Let $\pi(x)$ be the prime-counting function that gives the number of primes less than or equal to x, then
$\frac{1}{\log x+2}<\frac {\pi(x)}x<\frac{1}{\log x-4}.$

Check here for the proof.
By the Bayes formula, for a $m$ bit number $n$ ,let $N=2^m$ , we have

$\begin{equation} \begin{split} P(P|T_k)&=\frac{P(T_k|P)P(P)}{P(T_k|P)P(P)+P(T_k|C)P(C)}\\ &>\frac {\pi(N)/N}{\pi(N)/N+\frac1{4^k}(1-\pi(N)/N)}\\ &=\frac{1}{1+\frac1{4^N}(N/\pi(N)-1)}\\ &>\frac{1}{1+\frac1{4^k}(\log N+1)}\\ &>1-\frac{m}{4^{k-1}} \end{split} \end{equation}$

The last inequality require $k>\log m$ , 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 $1.8*10^{-24}$ . So one may choose $k\geq43\log \log n$ to guarantee the precision of the algorithm is suitable for practical use.

Remark 11.

It was show here that admitting the generalized Riemann hypothesis, if $a$ is a Miller-Rabin witness of $n$ , then $1\leq a\leq 2\log^2 n$ . This leads to a deterministic Miller-Rabin primality algorithm with polynomial time complexity $\tilde\sO(\log^4 n)$ ， 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 $\sO(\log s)$ times multiplication but the digit of $a$ increase to $\log^2 n$ digits, which has time complexity $\tilde\sO(\log^2 n)$ 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 $\tilde\sO(\log n)$ .

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

An introduction to Miller-Rabin algorithm and its Python

Miller-Rabin algorithm

1. Introduction

2. Math background for the Miller-Rabin primality test

1. Fermat’s primality test

2. Miller-Rabin primality test

3. Python implementation of Miller-Rabin primality test

CATALOG

FRIENDS