Skip to content

Instantly share code, notes, and snippets.

@tchaumeny
Created Jul 15, 2021
Embed
What would you like to do?
Implementation of Fermat and Miller-Rabin tests in Python, see https://lipsum.dev/2021-06-1-miller-rabin-openssl/
# Part of https://lipsum.dev/2021-06-1-miller-rabin-openssl/
# /!\ Warning /!\ — Not safe for real cryptographic usage
from collections import Counter
from random import randrange
import statistics
# Fermat
def fermat_test(n, max_witness=100):
witnesses = range(2, min(n, max_witness + 1))
for a in witnesses:
if pow(a, n - 1, n) != 1:
print(f"{n} is not a prime ({a} is a witness)")
return
print(f"{n} might be a prime")
fermat_test(4208312609) # 4208312609 is prime
fermat_test(4208312609**2)
fermat_test(5472940991761) # 5472940991761 = 199 * 241 * 863 * 132233
# Miller-Rabin
COMPOSITE = 'COMPOSITE'
NO_WITNESS = 'NO_WITNESS'
def miller_rabin_test(n, witnesses):
assert n >= 3
# Find greatest r s.t. 2^r | n - 1
r = 1
while (n - 1) % 2**(r + 1) == 0:
r += 1
d = (n - 1) // 2**r
for a in witnesses:
# Check first equality
x = pow(a, d, n)
if x in (1, n - 1):
continue
# Check successive squares
for i in range(1, r):
x = pow(x, 2, n)
if x == n - 1:
break
if x == 1:
return COMPOSITE
else:
return COMPOSITE
return NO_WITNESS
miller_rabin_test(4208312609**2, [2])
miller_rabin_test(5472940991761, [2])
miller_rabin_test(47197, [3])
def is_probable_prime(n, rounds=10):
if n == 1:
return False
if n in (2, 3):
return True
return miller_rabin_test(
n,
(randrange(2, n - 2) for _ in range(rounds))
) != COMPOSITE
is_probable_prime(4208312609)
is_probable_prime(4208312609**2)
# Primes generation stats
def gen_prime_trials(bits):
i = 1
while True:
rnd = randrange(2**(bits - 1) + 1, 2**bits, 2)
if is_probable_prime(rnd):
return i
i += 1
print("b\tTrials\tTrials / b")
for b in range(2, 64, 5):
avg = statistics.mean(gen_prime_trials(b) for _ in range(10**4))
print(f"{b}\t{float(avg):.3}\t{avg / b:.3}")
# Prime generation
SMALL_PRIMES = [
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149,
151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
479, 487, 491, 499
]
def has_small_factor(n):
for p in SMALL_PRIMES:
if p**2 > n:
return False
if n % p == 0:
return True
return False
def gen_prime(bits):
while True:
rnd = randrange(2**(bits - 1) + 1, 2**bits, 2)
delta = 0
while rnd + delta < 2**bits:
if has_small_factor(rnd + delta):
delta += 2
elif is_probable_prime(rnd + delta):
return rnd + delta
else:
break
gen_prime(2048)
c = Counter(gen_prime(10) for _ in range(10**4))
c.most_common(3)
c.most_common()[-3:]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment