from imp import reload
from math import sqrt,ceil
from random import randrange

""" These algorithms are taken from
    'Elementary Number Theory: Primes, Congruences, and Secrets'
    by Stein (with some modifications here and there).
"""

# Algorithm 1.1.13
def gcd(a,b):
    a,b = abs(a),abs(b)
    while b != 0:

        # a = q(b)+(r)
        q,r = divmod(a,b)

        # print('{} = {}({})+({})'.format(a,q,b,r))

        a,b = b,r

    return a

# Algorithm 2.3.7 (modified)
def xgcd(a,b):
    assert a>=0 and b>=0
    x,r = 1,0
    y,s = 0,1

    # print("{}=a({})+b({})".format(a,x,y))
    while b>0:
        q = a//b
        a,b = b,a-q*b
        x,r = r,x-q*r
        y,s = s,y-q*s

        # print("q={}\t{}=a({})+b({})".format(q,a,x,y))

    return a,x,y

# Algorithm 1.2.3
def prime_sieve(n,verbose=True):
    P = []
    X = range(2,n+1)

    while True:

        p=X[0]
        if p <= sqrt(n):
            P.append(p)
            X=[m for m in X if m%p != 0]
            if verbose:
                print("P = {}\nX = {}\n".format(P,X))

        else:
            P.extend(X)
            if verbose:
                print("P = {}\n".format(P))

            return P

def is_prime_miller_rabin(n,a=None,verbose=False):
    m=n-1
    k=0
    while m%2 == 0:
        m //= 2
        k += 1

    if not a:
        a = randrange(2,n)
    b = pow(a,m,n)
    if verbose:
        print("({0})^({1}) = {2} (mod {3})".format(a,m,b,n))
    if b == 1 or b == n-1:
        return True

    for r in range(1,k):
        c = pow(b,2,n)
        if verbose:
            print("({0})^({1}) = {2} (mod {3})".format(b,2,c,n))
        b = pow(b,2,n)
        if b == n-1:
            return True

    return False

def is_probable_prime(n,iterations=10):
    """ If k=iterations and this returns True the chance
        that 'n' is composite is less than 1/4**k.
    """
    for k in range(iterations):
        is_prime = is_prime_miller_rabin(n)
        if not is_prime:
            return False
    return True

def next_nice_prime(p):
    """ returns next prime Q larger that p such that
        (Q-1)/2 is also prime.
    """
    is_nice_prime = False
    while not is_nice_prime:
        p = next_probable_prime(p)
        q = (p-1)//2
        is_nice_prime = is_probable_prime(q)
    return p

def is_prime_lucas_lehmer(p):
    n = 2**p - 1
    s = 4
    for k in range(3,p+1):
        s = pow(s,2,n) - 2
    return s==0

def next_probable_prime(n,iterations=10):
    if n%2 == 0:
        n -= 1
    is_prime = False
    while not is_prime:
        n += 2
        is_prime = is_probable_prime(n,iterations=iterations)
    return n

def next_mersenne_prime(p):
    is_mersenne_prime = False
    while not is_mersenne_prime:
        p = next_probable_prime(p)
        is_mersenne_prime = is_prime_lucas_lehmer(p)
    return p

# p is prime but (p-1)//2 is not prime
p = 93450983094850938450983409611

# q is prime and (q-1)//2 is also prime
q = 93450983094850938450983409623


# RSA stuff

def encrypt(m,e,n):
    return pow(m,e,n)

def decrypt(c,d,n):
    return pow(c,d,n)

def rsa(bits,e=65537):
    """ Generates probable primes 'p' and 'q' then forms the
        product n=pq that will have binary length of 'bits'.
        Next, if 'e' hasn't already been chosen, a random 'e'
        relatively prime to phi(n)=(p-1)(q-1) is generated and
        'd' with ed=1 (mod phi(n)) is found. (e,d,n) are returned.
    """
    bound = 2**(bits//2 + 1)
    p=randrange(1,bound)
    p = next_probable_prime(p)

    q=randrange(1,bound)
    q = next_probable_prime(q)
    if p==q:
        q = next_probable_prime(q)

    n = p*q
    phi_n = (p-1)*(q-1)
    while True:
        # standard values for e are 3, 17, 65537
        if not e:
            e = randrange(2,phi_n)
        g,d,y = xgcd(e,phi_n)
        if g==1:
            while d<0:
                d += phi_n
            return e,d,n
        else:
            e = randrange(2,phi_n)

def encode(s):
    s = str(s)          # make input a string if it wasn't already
    L = [ord(c)*256**k for k,c in enumerate(s)]
    return sum(L)

def decode(n):
    n = int(n)          # make input an integer if it wasn't already
    v = []              # make a list for the characters
    while n != 0:
        n,r = divmod(n, 256)
        v.append(chr(r))
    return ''.join(v)

# Cracking RSA

# Cracking RSA if phi(n) is known
def crack_given_phi(n, phi_n):
    """ This factors n=pq when phi(n) is known.
        Note that phi(n)=(p-1)(q-1)=pq-(p+q)+1 so
        we know p+q = n-phi(n)+1. This means we
        have the coefficients of the polynomial
        (x-p)(x-q)=x^2-(p+q)x+pq and so can solve
        it with the quadratic formula.
    """
    #a= 1
    b = n - phi_n + 1           # -b in standard formula
    c = n
    d = int(sqrt(b**2 - 4*c))
    p = (b+d)//2
    q = (b-d)//2
    return p,q

# Factoring n=pq when p and q are 'close'
def crack_when_pq_close(n):
    """ Assume n=pq with p>q and define
                s=(p-q)/2,  t=(p+q)/2.
        Then
                n=t**2-s**2
        where s is small. Observe that
                t**2-n = s**2
        is a perfect square. We start with the guess that
                t = ceil(sqrt(n))
        and check to see if
                t**2 - n
        is a perfect square. If it isn't we bump up t by 1
        and check again. Eventually we'll find s and t.
        Then
                n = s+t
    """
    t = ceil(sqrt(n))
    while True:
        try:
            s2 = t**2-n
            s  = int(sqrt(s2))
            if s**2 == s2:
                p = t+s
                q = t-s
                return p,q
        except ValueError:
            pass
        t += 1

# Factoring n=pq when d is given
def crack_given_decrypt(n,e,d):
    m = e*d-1
    while True:
        if m%2 != 0:
            break
        divide_out = True
        for k in range(5):
            a = randrange(2,n)
            g = gcd(a,n)

            if g == 1:
                if pow(a,m//2,n) != 1:
                    divide_out = False

            # Or 1<g<n
            else:
                return g

        if divide_out:
            m = m//2
        else:
            break

    while True:
        a = randrange(2,n)
        g = gcd(pow(a,m//2,n)-1,n)
        if g != 1 and g != n:
            return g

# Example for 'crack_given_decrypt'
e = 22601762315966221465875845336488389513
d = 31940292321834506197902778067109010093
n = 268494924039590992469444675130990465673