RSA Algorithm to establish a random secret key. Go to cocalc to run the code. You can also install sage on your local machine.

reset()
        
# Return x^a mod n
def pow_mod(x,a,n) :
    L = []
    m = a
    while (m != 0) :
        L.append(m%2)
        m = m >> 1
    M = []
    for i in [0,..,len(L)-1] :
        if (i==0) : M.append(x)
        else : M.append( (M[-1]*M[-1]) % n )
    result = 1
    for i in [0,..,len(L)-1] :
        if ( L[i] == 1) :
            result *= M[i]
            result = result % n
    return result
        
        
# Return the GCD of integers a and b
def GCD(a,b) :
    while (b!=0) :
        r = (a%b)
        a=b
        b=r
    return a
        
        
# Return (gcd, alpha, beta) such that gcd = alpha.a + beta.b
def EGCD(a,b) :
    stack = []
    while (b!=0) :
        r = (a%b)
        m = (a-r)/b
        stack.append([m,-1,-1])
        a=b
        b=r
    gcd=a
    stack.append([-1,1,0])
    for i in range(len(stack)-2,-1,-1) :
        stack[i][1] = stack[i+1][2]
        stack[i][2] = stack[i+1][1] - stack[i][0]*stack[i+1][2]
    return (gcd, stack[0][1], stack[0][2])
        
        
# Basic Miller-Rabin with confidence 1-2^{-t}
def miller_rabin(N,t) :
    for i in [1,..,t] :
        x = 1 + ZZ.random_element(N-1)
        val = pow(x,N-1,N)
        #print val
        if ( val != 1 ) : return False
    return True
        
        
# Generate an n-bit prime
def gen_prime(n,t) :
    for i in [1,..,n*t] :
        p = 2^(n-1) + ZZ.random_element(2^(n-1))
        if miller_rabin(p,t) : return p
    return -1
        
        
# Find a random exponenet e
def find_exponent(phi_N, num) :
    for i in [1,..,num] :
        e = 2 + ZZ.random_element(phi_N-2)
        g = GCD(e,phi_N)
        if (g==1) : return e
    return -1
        
        
        
### MAIN CODE
### The RSA key-agreement protocol
        
# The number of bits in the primes
n=4096
        
# Generate the first prime
p = gen_prime(n,100)
print "p =", p
        
# Generate the second prime (different from the first one)
q = gen_prime(n,100)
while (q == p) : q = gen_prime(n,100)
print "q =", q
        
# Compute the modulus (the product of the two primes)
N = p*q
print "N =" , N
        
# Compute the size of the group (Z_N^*,X), i.e. \phi(N)
phi_N = (p-1) * (q-1)
print "phi(N) =", phi_N
        
# Compute a random exponent that is relatively prime to \phi(N)
e = find_exponent(phi_N,100)
print "e =", e

#Generate a random mask r
r = ZZ.random_element(N)
while (GCD(r,N) != 1) : r = ZZ.random_element(N)
print "r =", r
        
# Encode the mask r into y
y=pow_mod(r,e,N)
print "y =", y
        
# Calcualte the trapdoor d such that ed = 1 mod \phi(N)
d = EGCD(e,phi_N)[1]%phi_N
print "d =", d
        
# Decrypt y to recover r using the trapdoor d
r_tilde = pow_mod(y,d,N)
print "r_tilde =", r_tilde
        
# Check if the random mask was successfully recovered
print "Was Decoding Successful?" , r==r_tilde
        
# Compute the inverse of the decrypted result in preparation for RSA Decryption algorithm
inv_r_tilde = EGCD(r_tilde,N)[1]%N
print "Inverse of r_tilde =", inv_r_tilde