For cryptography, we often need large prime numbers. But it takes a lot of computation to check if a number is prime by trying to factor it.

Instead of trying to factor a number, we can use a Monte Carlo method to quickly check if a number is prime. The Miller-Rabin primality test is a probabilistic algorithm that can quickly check if a number is prime.

Algorithm

In [2]:
def miller_rabin(n: int, trials: int) -> bool:
    if n & 1 == 0: return False
    if n < 2: return False
    if n in (2, 3): return True

    s = 0
    d = n - 1

    while d & 1 == 0:
        s += 1
        d >>= 1

    for _ in range(trials):
        a = rng.randint(2, n - 2)
        x = pow(a, d, n)

        for _ in range(s):
            y = pow(x, 2, n)
            if y == 1 and x != 1 and x != n - 1: return False
            x = y

        if x != 1: return False

    return True

Let’s test this by computing the first few Mersenne primes.

In [3]:
for i in range(600):
    p = 2 ** i - 1
    if miller_rabin(p, 10):
        print(f"2^{i} - 1")
Out:
2^2 - 1
2^3 - 1
2^5 - 1
2^7 - 1
2^13 - 1
2^17 - 1
2^19 - 1
2^31 - 1
2^61 - 1
2^89 - 1
2^107 - 1
2^127 - 1
2^521 - 1

Looks like it works! Feel free to compare with the list of Mersenne primes.

Now let’s try to generate some random primes.

Generating random primes

In [4]:
def n_bit_prime(n: int) -> int:
    low = 2 ** (n - 1)
    high = 2 ** n - 1
    while True:
        p = rng.randint(low, high)
        p |= 1  # make sure it's odd
        if miller_rabin(p, 10):
            return p

n_bit_prime(64)
Out [4]:
10245293735733408889

To see if it’s really prime, we can try to factor it. This can be done easily with the factor command on Linux.

In [5]:
run("factor 69420")
Out:
69420: 2 2 3 5 13 89

Now let’s try to factor our random prime.

In [6]:
run(f"factor {n_bit_prime(64)}")
Out:
15031861708546515673: 15031861708546515673

Looks like it’s really prime!

Generating safe primes

For some use cases, we might need “safe” primes. A safe prime is a prime number p such that (p - 1) / 2 is also prime.

You can read more about safe primes on the Wikipedia page for Safe and Sophie Germain primes.

In [7]:
def safe_prime(n: int) -> int:
    while True:
        p = n_bit_prime(n)
        q = p * 2 + 1
        if miller_rabin(q, 10):
            return q

result = safe_prime(64)
run(f"factor {result}")
Out:
25662923996755684019: 25662923996755684019

Since this is a safe prime, (p - 1) / 2 is supposed to be prime as well. Let’s check that.

In [8]:
run(f"factor {(result - 1) // 2}")
Out:
12831461998377842009: 12831461998377842009

It is indeed prime! We have successfully generated a safe prime.

Useful links