#!/usr/bin/python
import random

def generateKey(bits=64):
	assert 64<=bits
	assert bits%4==0
	length=bits//4
	pswd=random.randrange(2**(bits-1), 2**bits)
	p,q=generatePrimePair(bits)
	n=p*q
        phi_n = (p-1)*(q-1)
        e = 65537 # for e = 3, k = 2!
	while not gcd(e,n):
		e+=2
	d=modInverse(e,phi_n)

        # Theorem: there exists a k < e, so that e * d = k * phi_n + 1 (mod 2^(bits/4))
        left = (e*d)%(2**length)
        for k in range(0, e):
            right = (k*phi_n + 1)%(2**length)
            if left == right:
                print "found k: %s" % k
        # is this a surprising result? Does that help at all with attacking RSA?

	return e,d,n

def RabinMillerWitness(test,possible):
	a,b,n=long(test%possible),possible-1,possible
	if a==1:
		return False
	A=a
	t=1L
	while t<=b:
		t<<=1
	#t=2**k, and t>b
	t>>=2
	while t:
		A=pow(A,2,n)
		if t&b:
			A=(A*a)%n
		if A==1:
			return False
		t>>=1
	return True

smallprimes = (3,5,7,11,13,17,19,23,29,31,37,41,43,
				47,53,59,61,67,71,73,79,83,89,97)

def getPrime(b,seed):
	bits=int(b)
	assert 64<=bits
	k=bits<<1
	possible=seed|1
	good=0
        while possible%12 != 11:
            possible+=2
	while not good:
		possible+=12
		good=1
		for i in smallprimes:
			if possible%i==0:
				good=0
				break
                for i in xrange(k):
                    test=random.randrange(2,possible)|1
                    if RabinMillerWitness(test,possible):
                        good=0
                        break
        return possible

def egcd(a,b):
        u, u1 = 1, 0
        v, v1 = 0, 1
        while b:
                q = a // b
                u, u1 = u1, u - q * u1
                v, v1 = v1, v - q * v1
                a, b = b, a - q * b
        return u, v, a

def gcd(a,b):
        a,b=(b,a) if a<b else (a,b)
        while b:
                a,b=b,a%b
        return a

def modInverse(e,n):
        return egcd(e,n)[0]%n

def generatePrimePair(seed,bits=64):
        assert 64<=bits
        assert bits%4==0
        seed1=random.randrange(2**(bits-2), 2**(bits-1))
        seed2=random.randrange(2**(bits-2), 2**(bits-1))
	p=getPrime(bits,seed1)
	q=getPrime(bits,seed2)
	return p,q

if __name__=="__main__":
	print generateKey(bits=256)