Implementation of Fermat and Miller-Rabin tests in Python, see https://lipsum.dev/2021-06-1-miller-rabin-openssl/

miller_rabin.py

# 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:]