import sympy, sys, random, math

sys.setrecursionlimit(10000)

def random_relatively_prime(n):
    while (True):
        candidate = random.randrange(n)
        if math.gcd(candidate, n) == 1:
            return candidate

# https://en.wikibooks.org/wiki/Algorithm_Implementation/Mathematics/Extended_Euclidean_algorithm

def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

def modular_inverse(a, m):
    g, x, y = egcd(a, m)
    return x % m

n_bits = 1024
p = sympy.randprime(3, 2 ** n_bits)
q = sympy.randprime(3, 2 ** n_bits)

N = p*q
phi = (p-1) * (q-1)

e = random_relatively_prime(phi)

d = modular_inverse(e, phi)

assert(math.gcd(e, phi) == 1)
assert((e * d) % phi == 1)

public_key = (e, N)
private_key = (d, N)

print("public key is " + str(public_key))
print("private key was " + str(private_key))

plaintext = 123456789123456789
ciphertext = pow(plaintext, public_key[0], public_key[1])
assert(plaintext == pow(ciphertext, private_key[0], private_key[1]))

signedtext = pow(plaintext, private_key[0], private_key[1])
assert(plaintext == pow(signedtext, public_key[0], public_key[1]))