next-two-adicity-prime.py

#!/usr/bin/env -S uv run --script
# /// script
# dependencies = [
# ]
# ///
# <https://gist.github.com/recmo/af5935fc8a40c54050a691301eaacdaa>
import random
import time

def miller_rabin(candidate, repeats=128):
    # Each Miller-Rabin iteration has at most a 1/4 chance of a false positive.
    # To achieve a false positive probability lower than 2^-256, we require:
    #       (1/4)^k < 2^-256
    # which simplifies to: k > 128.
    # Thus, performing 128 iterations is sufficient.
    d = candidate - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for _ in range(repeats):
        a = random.randrange(2, candidate - 1)
        x = pow(a, d, candidate)
        if x == 1 or x == candidate - 1:
            continue
        for _ in range(s - 1):
            x = pow(x, 2, candidate)
            if x == candidate - 1:
                break
        else:
            return False
    return True

assert(not miller_rabin(57))
assert(miller_rabin(2**31-1))
assert(miller_rabin(2**64-2**32+1))

k = 20
l = 128

while True:
    start_time = time.time()
    q = random.randrange(2**l, 2**(l+1))
    while True:
        p = q * 2**k + 1
        if miller_rabin(p):
            break
        q += 1
    elapsed_time = time.time() - start_time
    print(p, elapsed_time)