In [9]:
get_ipython().ast_node_interactivity = 'all'
import os
import random
import math

def bits_of_int(n):
    bits = 0
    while n:
        n //= 2
        bits += 1
    return bits

def random_int(high):
    bits = bits_of_int(high)
    mask = int('1' * bits, 2)

    nbyte = bits // 8
    if bits % 8:
        nbyte += 1
    
    while True:
        n = int.from_bytes(os.urandom(nbyte), "big")
        n &= mask
        if n > high:
            continue
        return n

def random_range(low, high):
    while True:
        n = random_int(high)
        if n >= low:
            return n
In [2]:
def random_nbit(n):
    mask = int('1' * n, 2)
    nbyte = n // 8
    if n % 8:
        nbyte += 1
    num = int.from_bytes(os.urandom(nbyte), "big")
    num &= mask
    
    # First bit always 1
    num |= 1 << (n - 1)
    
    # Always an odd number
    num |= 1
    
    return num
In [3]:
def miller_rabin(n, k):
    s, d = 0, n - 1
    
    while d % 2 == 0:
        d = d // 2
        s += 1
    
    for _ in range(k):
        a = random_range(2, n - 1)
        x = pow(a, d, n)
        
        if x == 1 or x == n - 1:
            continue
        
        for _ in range(s - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True
In [4]:
def random_prime_nbit(N):
    while True:
        n = random_nbit(N)
        if miller_rabin(n, 64):
            return n
In [5]:
def random_safe_prime_nbit(N):
    while True:
        p = random_prime_nbit(N)
        sp = 2 * p + 1
        if miller_rabin(sp, 64):
            return sp
In [6]:
random_safe_prime_nbit(64)
Out [6]:
31229329104201141107
In [7]:
def random_elgamal_pg(N):
    p = random_safe_prime_nbit(N)
    while True:
        g = pow(random_range(2, p - 2), 2, p)
        
        if g == 1:
            continue
        
        if g == 2:
            continue
            
        if (p - 1) % g == 0:
            continue
        
        ginv = pow(g, -1, p)
        if (p - 1) % ginv == 0:
            continue
        
        return p, g
In [10]:
def random_elgamal_key_pair(p, g):
    b = bits_of_int(p)
    priv = random_range(1 << (b - 1), p - 2)
    pub = pow(g, priv, p)
    return priv, pub
In [8]:
p, g = random_elgamal_pg(64)
In [11]:
priv, pub = random_elgamal_key_pair(p, g)