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

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.

for i in range(600): p = 2 ** i - 1 if miller_rabin(p, 10): print(f"2^{i} - 1")

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

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)

10030590156072197917

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

command on Linux.

run("factor 69420")

69420: 2 2 3 5 13 89

Now let’s try to factor our random prime.

run(f"factor {n_bit_prime(64)}")

16317264798946290947: 16317264798946290947

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.

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}")

19429457372028839327: 19429457372028839327

Since this is a safe prime, `(p - 1) / 2`

is supposed to be prime as well. Let’s check that.

run(f"factor {(result - 1) // 2}")

9714728686014419663: 9714728686014419663

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