##########################################################################
## File: proj2.py
## Author: Christopher S. Marron
## Date: 07 November 2014
##
## Utility functions for use in CMSC 426/626 Project 2
##########################################################################

import random

# Extended Euclidean Algorithm from Wikibooks.  Finds the
# gcd of the integer inputs.
#
#   a, b - integer inputs
#
# Returns b, x, y such that b is the gcd and b = a*x + b*y
#
def egcd(a, b):
    """
    Compute the GCD of a and b using the Extended Euclidean Algorithm.

    Returns [d, x, y] where d is the GCD of a and b and x and y are 
    two integers such that d = ax + by
    """

    x,y, u,v = 0,1, 1,0
    while a != 0:
        q, r = b//a, b%a
        m, n = x-u*q, y-v*q
        b,a, x,y, u,v = a,r, u,v, m,n
    return b, x, y

# Modular inverse from Wikibooks.  Given a and m, determines
# the inverse of a mod m.
#
#   a - integer to find the inverse of
#   m - modulus
#
def modinv(a, m):
    """
    Modular inverse computation.

    Returns the inverse of a modulo m. Computes the inverse
    using the Extended Euclidean Algorithm.
    """

    g, x, y = egcd(a, m)
    if g != 1:
        return None  # modular inverse does not exist
    else:
        return x % m

 
_mrpt_num_trials = 5 # number of bases to test for Miller-Rabin

# Miller-Rabin primality test from rosettacode.org
#
#    n - integer to test
#
# Returns True if n is a probable prime
#
def is_probable_prime(n):
    """
    Miller-Rabin primality test.
 
    A return value of False means n is certainly not prime. A return value of
    True means n is very likely a prime.
    """
    assert n >= 2
    # special case 2
    if n == 2:
        return True
    # ensure n is odd
    if n % 2 == 0:
        return False
    # write n-1 as 2**s * d
    # repeatedly try to divide n-1 by 2
    s = 0
    d = n-1
    while True:
        quotient, remainder = divmod(d, 2)
        if remainder == 1:
            break
        s += 1
        d = quotient
    assert(2**s * d == n-1)
 
    # test the base a to see whether it is a witness for the compositeness of n
    def try_composite(a):
        if pow(a, d, n) == 1:
            return False
        for i in range(s):
            if pow(a, 2**i * d, n) == n-1:
                return False
        return True # n is definitely composite
 
    for i in range(_mrpt_num_trials):
        a = random.randrange(2, n)
        if try_composite(a):
            return False
 
    return True # no base tested showed n as composite