from math import sqrt,log

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

def gcd(a,b):
    a,b=abs(a),abs(b)
    while b!=0:
        a,b = b,a%b
    return a

##################################################
## Linear Equations Modulo $n$
##################################################

def sign(n):
    """
    Returns -1 if n<0 and +1 otherwise.
    """
    if n<0:
        return -1
    return 1

def xgcd(a, b):
    """
    Returns g, x, y such that g = x*a + y*b = gcd(a,b).
    Input:
        a -- an integer
        b -- an integer
    Output:
        g -- an integer, the gcd of a and b
        x -- an integer
        y -- an integer
    """
    a,x_sign = abs(a),sign(a)
    b,y_sign = abs(b),sign(b)

    if a == 0:
        return (b,0,y_sign)
    if b == 0:
        return (a,x_sign,0)

    x = 1; y = 0; r = 0; s = 1
    while b != 0:
        q = a//b
        (a, x, y, b, r, s) = (b, r, s, a-q*b, x-q*r, y-q*s)
    return (a, x*x_sign, y*y_sign)

def inversemod(a, n):
    """
    Returns the inverse of a modulo n, normalized to
    lie between 0 and n-1.  If a is not coprime to n,
    raise an exception (this will be useful later for 
    the elliptic curve factorization method).
    Input:
        a -- an integer coprime to n
        n -- a positive integer
    Output:
        an integer between 0 and n-1.
    """
    g, x, y = xgcd(a, n)
    if g != 1:
        raise ZeroDivisionError((a,n))
    assert g == 1, "a must be coprime to n."
    return x%n


##################################################
## Elliptic Curve Addition
##################################################

def ellcurve_add(E, P1, P2):
    """
    Returns the sum of P1 and P2 on the elliptic 
    curve E.
    Input:
         E -- an elliptic curve over Z/pZ, given by a 
              triple of integers (a, b, p), with p odd.
         P1 --a pair of integers (x, y) or the 
              string "Identity".
         P2 -- same type as P1
    Output:
         P3 -- same type as P1
    """
    a, b, p = E
    assert p > 2, "p must be odd."

    if P1 == "Identity":
        P3 = P2

    elif P2 == "Identity":
        P3 = P1

    else:
        x1, y1 = P1
        x2, y2 = P2
        x1 %= p; y1 %= p; x2 %= p; y2 %= p

        if x1 == x2:
            if (y1+y2)%p == 0:          # P1 ==-P2
                return "Identity"

            else:                       # P1 == P2
                lam = (3*x1**2+a) * inversemod(2*y1,p)

        else:                           # P1 != P2
            lam = (y1-y2) * inversemod(x1-x2, p)

        x3 = lam**2 - x1 - x2
        y3 = -(y1 + lam*(x3-x1))
        P3 = (x3%p, y3%p)

    return P3