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)
18332724092051638669
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)}")
13978716916633231963: 13978716916633231963
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}")
19936227763520711003: 19936227763520711003
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}")
9968113881760355501: 9968113881760355501
It is indeed prime! We have successfully generated a safe prime.