def double_affine(P, a):
    """
    INPUT:
    - `P`: point in affine coordinates
    - `a`: curve coefficient in y^2 = x^3 + a + b

   Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw.html
    """
    x1, y1 = P[0], P[1]
    lambd = (3*x1**2+a)/(2*y1)
    x3 = lambd**2-x1-x1
    y3 = lambd*(x1-x3)-y1
    return P.curve()((x3,y3,1))

def add_affine(P, Q):
    """
    INPUT:
    - `P`: point in affine coordinates
    - `Q`: point in affine coordinates
    - `a`: curve coefficient in y^2 = x^3 + a + b

   Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw.html
    """
    x1, y1 = P[0], P[1]
    x2, y2 = Q[0], Q[1]
    lambd = (y2-y1)/(x2-x1)
    lambd2 = lambd**2
    x3 = lambd2-x1-x2
    y3 = lambd*(x1-x3)-y1
    return P.curve()((x3,y3,1))

def group_law_affine(P, Q):
    E = P.curve()
    assert E.a2() == 0 and E.a1() == 0 and E.a3() == 0 # short Weierstrass form
    x1, y1 = P[0], P[1]
    x2, y2 = Q[0], Q[1]
    if x1 == x2 and y1 == -y2: # P == -Q, P+Q = P-P = 0
        return P.curve()(0) # the point at infinity
    if x1 == x2 and y1 == y2:
        return double_affine(P, E.a4())
    else:
        return add_affine(P, Q)

def double_projective(P, a):
    """
    INPUT:
    - `P`: point in projective coordinates
    - `a`: curve coefficient in Y^2*Z = X^3 + a*X*Z^2 + b*Z^3

    Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#doubling-dbl-1998-cmo-2
    """
    X1,Y1,Z1 = P # gets the projective coordinates of P
    w = a*Z1**2+3*X1**2
    s = Y1*Z1
    ss = s**2
    sss = s*ss
    R = Y1*s
    B = X1*R
    h = w**2-8*B
    X3 = 2*h*s
    Y3 = w*(4*B-h)-8*R**2
    Z3 = 8*sss
    # at the end, cast the point to the curve
    return P.curve()((X3, Y3, Z3))

def add_projective(P, Q):
    """
    Addition in projective coordinates on a curve
    in reduced Weierstrass form Y^2*Z = X^3 + a*X*Z^2 + b*Z^3
    INPUT:
    - `P`: point (X1,Y1,Z1) in projective coordinates
    - `Q`: point (X2,Y2,Z2) in projective coordinates

    Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-1998-cmo-2
    """
    X1,Y1,Z1 = P # gets the projective coordinates of P
    X2,Y2,Z2 = Q # gets the projective coordinates of Q
    Y1Z2 = Y1*Z2
    X1Z2 = X1*Z2
    Z1Z2 = Z1*Z2
    u = Y2*Z1-Y1Z2
    uu = u**2
    v = X2*Z1-X1Z2
    vv = v**2
    vvv = v*vv
    R = vv*X1Z2
    A = uu*Z1Z2-vvv-2*R
    X3 = v*A
    Y3 = u*(R-A)-vvv*Y1Z2
    Z3 = vvv*Z1Z2
    return P.curve()(X3, Y3, Z3)

def double_jacobian(P, a):
    """
    INPUT:
    - `P`: point in jacobian coordinates
    - `a`: curve coefficient in Y^2 = X^3 + a*X*Z^4 + b*Z^6

    Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-1998-cmo-2
    It is not possible to cast at the end because Jacobian coordinates are not native in SageMath
    """
    X1,Y1,Z1 = P # gets the Jacobian coordinates of P
    XX = X1**2
    YY = Y1**2
    ZZ = Z1**2
    S = 4*X1*YY
    M = 3*XX+a*ZZ**2
    T = M**2-2*S
    X3 = T
    Y3 = M*(S-T)-8*YY**2
    Z3 = 2*Y1*Z1
    return (X3, Y3, Z3)

def add_jacobian(P, Q):
    """
    INPUT:
    - `P`: point (X1,Y1,Z1) in Jacobian coordinates
    - `Q`: point (X2,Y2,Z2) in Jacobian coordinates

    Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-1998-cmo-2
    It is not possible to cast at the end because Jacobian coordinates are not native in SageMath
    """
    X1,Y1,Z1 = P # gets the Jacobian coordinates of P
    X2,Y2,Z2 = Q # gets the Jacobian coordinates of Q
    Z1Z1 = Z1**2
    Z2Z2 = Z2**2
    U1 = X1*Z2Z2
    U2 = X2*Z1Z1
    S1 = Y1*Z2*Z2Z2
    S2 = Y2*Z1*Z1Z1
    H = U2-U1
    HH = H2
    HHH = H*HH
    r = S2-S1
    V = U1*HH
    X3 = r2-HHH-2*V
    Y3 = r*(V-X3)-S1*HHH
    Z3 = Z1*Z2*H
    return (X3, Y3, Z3)

def addition_affine_edwards(P, Q, c, d):
    """
    Edwards form u^2 + v^2 = c^2*(1 + d*u^2 + v^2)

    INPUT:
    - `P`: point (u1,v1) in affine Edwards coordinates
    - `Q`: point (u2,v2) in affine Edwards coordinates
    """

    u1, v1 = P[0], P[1]
    u2, v2 = Q[0], Q[1]

    u3 = (u1*v2 + u2*v1)/(c*(1+d*u1*u2*v1*v2))
    v3 = (v1*v2 - u1*u2)/(c*(1 - d*u1*u2*v1*v2))

    return (u3, v3)