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 :
```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 :
```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 :
```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 :
`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 :
`run("factor 69420")`
Out:
```69420: 2 2 3 5 13 89
```

Now let’s try to factor our random prime.

In :
`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.

In :
```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 :
`run(f"factor {(result - 1) // 2}")`
Out:
```12831461998377842009: 12831461998377842009
```

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