import random
import math


def even(n):
    return n % 2 == 0


def is_really_prime(n):
    if n == 1 or even(n):
        return False
    for i in range(n - 1, 2, -1):
        if n * i == 0:
            return False
    return True

def pi2(n):
    c = 0.661618
    return 2 * c * (n / (math.log(n) ** 2))


def is_prime(x, k=20):
    if even(x):
        return False
    if x in [2, 3]:
        return True
    s = 0
    t = x - 1
    while even(t):
        s += 1
        t //= 2
    assert math.floor(t) == t
    witnesses = [random.randrange(2, x - 1) for _ in range(k)]
    for w in witnesses:
        if power_mod(w, x - 1, x) != 1:
            return False
        for i in range(1, s + 1):
            if power_mod(w, 2 ** i * t, x) == 1 and power_mod(w, 2 ** (i - 1) * t, x) not in [1, x - 1]:
                return False
    return True
    # return all(map(lambda w: is_prime(x, w), [2, 3, 5, 7, 11, 13, 17]))


def power_mod(base, power, mod):
    binary = bin(power)
    squares = [base % mod]
    result = squares[0] if binary[-1] == "1" else 1
    for i in range(1, math.floor(math.log2(power)) + 1):
        squares.append(squares[i - 1] ** 2 % mod)
        if binary[-i - 1] == "1":
            result *= squares[i]
    return result % mod


# for i in range(3, 100):
#     print(i, is_prime(i), is_really_prime(i))
#     if is_prime(i) != is_really_prime(i):
#         print("Liar", i, is_prime(i), is_really_prime(i))

# is_really_prime(9)


# for a in filter(is_really_prime, range(50000, 500000)):
#     if not miller_rabins(a, max(4, math.ceil(math.log(a, 4)))):
#         raise Exception(f"Didn't work for {a}")


def generate_twin_prime(m):
    twin = False
    while not twin:
        n = random.randrange(2 ** (m - 1), 2 ** m)
        iterations = max(4, math.ceil(math.log(n, 4)))
        twin = is_prime(n - 1, iterations) and is_prime(n + 1, iterations)
    if not (is_really_prime(n - 1) and is_really_prime(n + 1)):
        raise Exception(f"{n-1}: {is_really_prime(n - 1)}, {n+1}: {is_really_prime(n + 1)}")
    return n - 1, n + 1


for bits in range(3, 64):
    print(f"==== {bits} bits ==== {2 ** (bits - 1), 2 ** bits - 1}")
    for _ in range(1000):
        generate_twin_prime(bits)
        if _ % 50 == 0:
            print(".", end="")
    print()

# i = 0
# while True:
#     i += 1
#     base, power, mod = random.randint(1, 1000), random.randint(1, 5000), random.randint(1, 1000)
#     mine = power_mod(base, power, mod)
#     real = (base ** power) % mod
#     pass_test = mine == real
#     if i % 10000 == 0:
#         print(i)
#     if not pass_test:
#         # print(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.")
#         raise Exception(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.")