from sage.all_cmdline import *   # import sage library

from sage.rings.integer_ring import ZZ
from sage.rings.finite_rings.finite_field_constructor import FiniteField, GF
from sage.schemes.elliptic_curves.constructor import EllipticCurve
from sage.arith.misc import GCD, gcd

def dbl_affine_montgomery(R, A, B):
    """
    Doubling on a Montgomery curve E: B*y^2 = x^3 + A*x^2 + x

    INPUT:
    - `R`: point in affine coordinates (x, y, 1) or (0, y, 0) (point at infinity)
    - `A`: curve coefficient, not 2 or -2
    - `B`: curve coefficient, non-zero

    RETURN: Q(x3, y3, 1) = 2*R in affine coordinates, or (0, 1, 0)
    PROBLEM: the point at infinity has no affine representation
    The third coordinate indicates whether the point is a infinity
    http://www.hyperelliptic.org/EFD/g1p/auto-montgom.html
    """
    Fq = A.parent()
    Fq1 = Fq(1)
    Fq0 = Fq(0)
    x1, y1, z1 = R[0], R[1], R[2]
    if z1 != Fq0 and z1 != Fq1:
        raise ValueError("Error dbl_affine_montgomery the input point R is not in affine coordinates (x1, y1, 1): R = {}".format(R))
    if y1 == 0:
        return (Fq0, Fq1, Fq0) # point at infinity
    inv = 1/(2 * B * y1)                   # non-zero
    l = (3 * x1**2 + 2 * A * x1 + 1) * inv # lambda
    l2 = l**2                              # lambda^2
    x3 = B * l2 - A - x1 - x1
    y3 = l * (x1 - x3) - y1
    return (x3, y3, Fq1)

def add_affine_montgomery(P, R, A, B):
    """
    Addition on a Montgomery curve E: B*y^2 = x^3 + A*x^2 + x

    INPUT:
    - `P`: point in affine coordinates (x1, y1, 1) or (0, y1, 0) (point at infinity)
    - `R`: point in affine coordinates (x2, y2, 1) or (0, y2, 0) (point at infinity)
    - `A`: curve coefficient, not 2 or -2
    - `B`: curve coefficient, non-zero

    RETURN: Q(x3, y3, 1) = R+P in affine coordinates, or (0, 1, 0)
    PROBLEM: the point at infinity has no affine representation
    The third coordinate indicates whether the point is a infinity
    http://www.hyperelliptic.org/EFD/g1p/auto-montgom.html
    """
    Fq = A.parent()
    Fq1 = Fq(1)
    Fq0 = Fq(0)
    x1, y1, z1 = P[0], P[1], P[2]
    x2, y2, z2 = R[0], R[1], R[2]
    if z1 != Fq0 and z1 != Fq1:
        raise ValueError("Error add_affine_montgomery the input point P is not in affine coordinates (x1, y1, 1) or at infinity (0, y1, 0): P = {}".format(P))
    if z2 != Fq0 and z2 != Fq1:
        raise ValueError("Error add_affine_montgomery the input point R is not in affine coordinates (x2, y2, 1) or at infinity (0, y2, 0): R = {}".format(R))
    if z1 == Fq0: # P is at infinity, return R
        return (x2, y2, z2)
    if z2 == Fq0: # R is at infinity, return P
        return (x1, y1, z1)
    # now the addition
    if (x1 == x2) and (y1 == -y2):# points are opposite
        return (Fq0, Fq1, Fq0) # point at infinity
    if (x1 == x2) and (y1 == y2): # points are the same: P+P = [2]P
        return dbl_affine_montgomery(P, A, B)
    # here we can assume that (y1 != 0 or y2 != 0) or (x1 != x2)
    inv = 1/(x2-x1)         # non-zero
    l = (y2-y1) * inv       # lambda
    l2 = l**2               # lambda^2
    x3 = B*l2 - A - x1 - x2
    y3 = l*(x1 - x3) - y1
    return (x3, y3, 1)

def dbl_proj_montgomery(R, A, B):
    """
    Doubling on a Montgomery curve E: B*y^2*Z = x^3 + A*x^2*Z + x*Z^2

    INPUT:
    - `R`: point in projective coordinates (X,Y,Z), if Z = 0 the point is at infinity
    - `A`: curve coefficient, not 2 or -2
    - `B`: curve coefficient, non-zero

    RETURN: Q(X3, Y3, Z3) = 2*R in projective coordinates
    """
    Fq = A.parent()
    Fq1 = Fq(1)
    Fq0 = Fq(0)
    X1, Y1, Z1 = R[0], R[1], R[2]
    if Z1 == Fq0 or Y1 == Fq0: # R is at infinity, or has order 2, return O
        return (Fq0, Fq1, Fq0)
    s = 2 * Y1 * Z1                  # non-zero
    t = X1 * Z1
    u = 3*X1**2 + 2*t * A + Z1**2
    uu = u**2
    sB = s * B
    v = s * sB # s**2 * B
    Y2 = 2 * Y1**2
    Y4 = 2 * Y2 # 4*Y1**2
    R = t * Y4 * B
    S = uu - A*v - 2*R
    X3 = S * sB
    Y3 = u * (R - S) - Y2 * sB**2
    Z3 = v * sB
    return (X3, Y3, Z3)

def mixed_add_proj_montgomery(R, P, A, B):
    """
    Mixed Addition on a Montgomery curve E: B*y^2*Z = x^3 + A*x^2*Z + x*Z^2

    INPUT:
    - `R`: point in projective coordinates (X1, Y1, Z1), if Z1 = 0 the point is at infinity
    - `P`: point in affine coordinates (x2, y2, 1) or (0, y2, 0)
    - `A`: curve coefficient, not 2 or -2
    - `B`: curve coefficient, non-zero

    RETURN: Q(X3, Y3, Z3) = R+P in projective coordinates
    """
    Fq = A.parent()
    Fq1 = Fq(1)
    Fq0 = Fq(0)
    X1, Y1, Z1 = R[0], R[1], R[2]
    x2, y2, z2 = P[0], P[1], P[2]
    if z2 != Fq0 and z2 != Fq1:
        raise ValueError("Error mixed_add_proj_montgomery the input point P is not in affine coordinates (x2, y2, 1) or at infinity (0, y2, 0): P = {}".format(P))
    if Z1 == Fq0: # R is at infinity, return P
        return P
    if z2 == Fq0: # P is at infinity, return R
        return R
    # now the addition
    u = y2 * Z1 - Y1
    v = x2 * Z1 - X1
    if (v == 0) and (u == -2*Y1): # x2 == X1/Z1 and Y1/Z1 == -y2 <=> -Y1 == y2*Z1 <=> y2*Z1 - Y1 == -2*Y1, points are opposite
        return (Fq0, Fq1, Fq0)
    if (v == 0) and (u == 0): # x2 == X1/Z1 and Y1/Z1 == y2 points are equal
        return dbl_proj_montgomery(P, A, B)
    uu = u**2
    vv = v**2
    vvv = v * vv
    R = vv * X1
    S = B * uu * Z1 - vvv - 2*R - A * Z1 * vv
    X3 = v * S
    Y3 = u * (R-S) - vvv * Y1
    Z3 = vvv * Z1
    return (X3, Y3, Z3)

def double_and_add_proj_montgomery(m, P, A, B):
    """
    scalar multiplication [m]*P in Montgomery curve B*y^2 = x^3 + A*x^2 + x

    INPUT:
    - `P` in affine coordinates (x1, y1, 1) or infinity (0, y1, 0)
    - `m` non-zero scalar

    RETURN: point R = [m]*P in projective coordinates
    """
    Fq = A.parent()
    Fq1 = Fq(1)
    Fq0 = Fq(0)
    Z1 = P[2]
    if m == 0 or Z1 == Fq0:
        return (Fq0, Fq1, Fq0)
    if Z1 != Fq0 and Z1 != Fq1:
        raise ValueError("Error double_and_add_proj_montgomery the input point P is not in affine coordinates (x2, y2, 1) or at infinity (0, y2, 0): P = {}".format(P))
    if m < 0:
        m = -m
        P = (P[0], -P[1], Fq1)
    m = ZZ(m)
    b = m.bits()
    n = len(b)
    R = (P[0], P[1], Fq1)
    for i in range(n - 2, -1, -1):
        R = dbl_proj_montgomery(R, A, B)
        if b[i] == 1:
            R = mixed_add_proj_montgomery(R, P, A, B)
    return R

def add_affine_x_only(x1, x2, x4):
    """
    x-only addition for the Montgomery ladder on E: B*y^2 = x^3 + A*x^2 + x

    INPUT:
    - `x1`: x-coordinate of point P in affine coordinates (x1, y1, 1)
    - `x2`: x-coordinate of point Q in affine coordinates (x2, y2, 1)
    - `x4`: x-coordinate of point (P-Q) in affine coordinates (x4, y4, 1)

    RETURN: x3 the x-coordinate of (P+Q)

    formula: x3*x4*(x1 - x2)^2 = (x1*x2 - 1)^2
    The function assumes: x1 != x2, x1 or x2 non-zero, and x4 is non-zero.
    It does not test for the consistency of the inputs.
    """
    x3 = (x1*x2 - 1)**2 / (x4 * (x1 - x2)**2)
    return x3

def double_affine_x_only(x1, A):
    """
    x-only doubling for the Montgomery ladder on E: B*y^2 = x^3 + A*x^2 + x

    INPUT:
    - `x1`: x-coordinate of point P in affine coordinates (x1, y1, 1)
    - `A`: curve coefficient: B*y^2 = x^3 + A*x^2 + x

    RETURN: x2 the x-coordinate of [2]P

    formula: 4*x1*x3*(x1**2 + A*x1 + 1) = (x1**2 - 1)**2
    """
    x3 = (x1**2 - 1)**2 / (4 * x1 * (x1**2 + A*x1 + 1))
    return x3

def add_proj_xz_only(X1, Z1, X2, Z2, X4, Z4):
    """
    Addition in Projective X,Z-only on a Montgomery curve E: B*Y^2*Z = X^3 + A*X^2*Z + X*Z^2

    INPUT:
    - `X1`, `Z1`: point P in projective (X,Z)-coordinates (X1, Z1), if Z1 = 0 the point is at infinity
    - `X2`, `Z2`: point R in projective (X,Z)-coordinates (X2, Z2), if Z2 = 0 the point is at infinity
    - `X4`, `Z4`: point P-R in projective (X,Z)-coordinates (X4, Z4), if Z4 = 0 the point is at infinity

    RETURN: (X3, Z3) = R+P in projective (X,Z)-coordinates
    Source: http://www.hyperelliptic.org/EFD/g1p/auto-montgom-xz.html#diffadd-dadd-1987-m
    X5 = Z1*(X2*X3-Z2*Z3)^2
    Z5 = X1*(X2*Z3-Z2*X3)^2
    """
    if Z1 == 0 and Z2 == 0:
        return (0, 0) # corresponds to (0, 1, 0) in full (X, Y, Z) coordinates
    if Z1 == 0:
        return (X2, Z2)
    if Z2 == 0:
        return (X1, Z1)
    T = X1 * Z2 - X2 * Z1
    # we just assume that T is non-zero
    X3 = Z4 * (X1 * X2 - Z1 * Z2)**2
    Z3 = X4 * T**2
    return X3, Z3

def mixed_add_proj_xz_only(X1, Z1, X2, Z2, x4):
    """
    Mixed Addition in projective (X,Z)-coordinates on a Montgomery curve E: B*Y^2*Z = X^3 + A*X^2*Z + X*Z^2

    INPUT:
    - `X1`, `Z1`: point P in projective (X,Z)-coordinates (X1, Z1), if Z1 = 0 the point is at infinity
    - `X2`, `Z2`: point R in projective (X,Z)-coordinates (X2, Z2), if Z2 = 0 the point is at infinity
    - `x4`: point P-R in affine x-only coordinates (not at infinity)

    RETURN: (X3, Z3) = R+P in projective (X,Z)-coordinates
    Source: http://www.hyperelliptic.org/EFD/g1p/auto-montgom-xz.html#diffadd-dadd-1987-m with Z4 = 1

    Exceptional cases:
    if X1/Z1 = X2/Z2 and Z1*Z2 != 0 <=> X1*Z2 - X2*Z1 = 0 then it should raise and exception
    X3 = (X1 * X2 - Z1 * Z2)**2
    Z3 = x4 * (X1 * Z2 - X2 * Z1)**2
    then Z3 contains (X1 * Z2 - X2 * Z1) to be tested
    """
    if Z1 == 0 and Z2 == 0:
        return (0, 0) # corresponds to (0, 1, 0) in full (X, Y, Z) coordinates
    if Z1 == 0:
        return (X2, Z2)
    if Z2 == 0:
        return (X1, Z1)
    T = X1 * Z2 - X2 * Z1
    # we just assume that T is non-zero
    X3 = (X1 * X2 - Z1 * Z2)**2
    Z3 = x4 * T**2
    return X3, Z3

def double_proj_xz_only(X1, Z1, A):
    """
    Doubling in projective (X,Z)-coordinates on a Montgomery curve E: B*Y^2*Z = X^3 + A*X^2*Z + X*Z^2

    INPUT:
    - `X1`, `Z1`: point P in projective (X,Z)-coordinates (X1, Z1), if Z1 = 0 the point is at infinity
    - `A`: curve coefficient: B*y^2 = x^3 + A*x^2 + x

    RETURN: (X2,Z2) the (X,Z)-coordinates of [2]P

    x3 = (x1**2 - 1)**2/(4*x1*(x1**2 + A*x1 + 1))
    X3/Z3 = Z1 * (X1**2 - Z1**2)**2/(4*X1 * (X1**2 + A*X1 * Z1 + Z1**2))
    if the point is a 2-torsion point (0,0,1) or (a,0,1) it will return (0,0) in (X,Z)
    """
    if Z1 == 0 or X1 == 0:
        return (0, 0)
    a = (X1**2 + A*X1 * Z1 + Z1**2)
    if a == 0: # 2-torsion point
        return (0, 0)
    X3 = (X1**2 - Z1**2)**2
    Z3 = 4*X1 * Z1 * a
    return X3, Z3

# Montgomery ladder in affine and projective coordinates
def montgomery_ladder_affine_x_only(m, P, A):
    """
    Montgomery ladder in affine coordinates

    INPUT:
    - `m`: scalar (integer)
    - `P`: point in projective coordinates (X, Y, Z)
    - `A`: curve coefficient
    """
    if m == 0:
        return (0, 0)
    if m < 0:
        m = -m
    mbits = ZZ(m).bits()
    xP = P[0]
    x0, x1 = xP, double_affine_x_only(xP, A)
    for i in range(len(mbits)-2, -1, -1):
        bi = mbits[i]
        if bi == 0:
            (x0, x1) = double_affine_x_only(x0, A), add_affine_x_only(x0, x1, xP)
        else:
            (x0, x1) = add_affine_x_only(x0, x1, xP), double_affine_x_only(x1, A)
    return x0
            
def montgomery_ladder_proj_xz_only(m, P, A):
    """
    Montgomery ladder in projective coordinates

    INPUT:
    - `m`: scalar (integer)
    - `P`: point in projective coordinates (X, Y, 1) with Z=1
    - `A`: curve coefficient
    """
    if m == 0:
        return (0, 0)
    if m < 0:
        m = -m
    mbits = ZZ(m).bits()
    XP = P[0]
    ZP = P[2]
    assert ZP == 1
    X0 = (XP, ZP)
    X1 = double_proj_xz_only(XP, ZP, A)
    for i in range(len(mbits)-2, -1, -1):
        bi = mbits[i]
        if bi == 0:
            (X0, X1) = double_proj_xz_only(X0[0], X0[1], A), add_proj_xz_only(X0[0], X0[1], X1[0], X1[1], XP, ZP)
        else:
            (X0, X1) = add_proj_xz_only(X0[0], X0[1], X1[0], X1[1], XP, ZP), double_proj_xz_only(X1[0], X1[1], A)
    return X0

################################################################################
## TESTS FOR THE HAND-IN
################################################################################

def test_add_affine_x_only(E, A, B):
    """
    test affine addition of points in Montgomery form and Montgomery x-only coordinates
    E is in SageMath-compatible Montgomery form: y^2 = x^3 + A/B*x^2 + 1/B^2*x
    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    It is assumed that the inputs are consistent, that is the points are distincts, non-opposite, and non-zero
    """
    # check that E is in SageMath compatible Montgomery form
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    # first, test addition of a torsion point with another distinct random point
    list_points = get_torsion_points(E, A, B)
    Fq = E.base_field()
    Fq1 = Fq(1)
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        Sw = E.random_element()
        while Sw[0] == R[0]:
            Sw = E.random_element()
        # now R and S are consistent inputs
        Rw = E(R[0]/B, R[1]/B, R[2])  # long Weierstrass coordinates
        S = (Sw[0]*B, Sw[1]*B, Sw[2]) # Montgomery coordinates
        Tw = Rw - Sw                  # long Weierstrass coordinates
        T = (Tw[0]*B, Tw[1]*B, Tw[2])
        x4 = T[0]
        x3 = add_affine_x_only(R[0], S[0], x4)
        Qw = Rw + Sw
        Q = (Qw[0]*B, Qw[1]*B, Qw[2])
        ok = Q[0] == x3
        if not ok:
            print("Error add_affine_x_only")
            print("input: R={}".format(R))
            print("input: S={}".format(S))
            print("input: R-S={}".format(T))
            print("exp out: x3={} (and y3={}), z3={}".format(Q[0], Q[1], Q[2]))
            print("got out: x3={}".format(x3))
        i = i+1
    i = 0
    while ok and i < 10:
        distinct=False
        while not distinct:
            Rw = E.random_element() # long Weierstrass coordinates
            Sw = E.random_element() # long Weierstrass coordinates
            distinct = Rw[0] != Sw[0]
        Tw = Rw - Sw
        R = (Rw[0]*B, Rw[1]*B, Fq1) # Montgomery coordinates
        S = (Sw[0]*B, Sw[1]*B, Fq1) # Montgomery coordinates
        T = (Tw[0]*B, Tw[1]*B, Fq1) # Montgomery coordinates
        x3 = add_affine_x_only(R[0], S[0], T[0]) # Montgomery coordinates
        Qw = Rw+Sw # long Weierstrass coordinates
        Q = (Qw[0]*B, Qw[1]*B, Qw[2])
        ok = x3 == Q[0]
        if not ok:
            print("Error add_affine_x_only")
            print("input: x1={}".format(R[0]))
            print("input: x2={}".format(S[0]))
            print("input: x4={}".format(T[0]))
            print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
            print("got out: x3={}".format(x3))
        i = i+1
    print("test add_affine_x_only with B={}: {}".format(B, ok))
    return ok

def test_double_affine_x_only(E, A, B):
    """
    test affine doubling of points in Montgomery form and Montgomery x-only coordinates
    E is in SageMath-compatible Montgomery form: y^2 = x^3 + A/B*x^2 + 1/B^2*x
    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    list_points = get_torsion_points(E, A, B)
    Fq = E.base_field()
    Fq1 = Fq(1)
    # first test with the torsion points, but not the 2-torsion points
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        # the double_affine_x_only function assumes that x1 is non-zero, and (x1^2 + A*x1 + 1) is non-zero
        if R[1] == 0:
            i=i+1
            continue # it will skip the rest of the iteration and directly go to the next i
        Rw = E(R[0]/B, R[1]/B) # long Weierstrass coordinates
        Qw = 2*Rw # long Weierstrass coordinates
        Q = (Qw[0]*B, Qw[1]*B, Qw[2])
        x3 = double_affine_x_only(R[0], A)
        ok = Q[0] == x3
        if not ok:
            print("Error double_affine_x_only")
            print("input: R[0]={}".format(R[0]))
            print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
            print("got out: x3={}".format(x3))
        i = i+1
    i = 0
    while ok and i < 10:
        order2=True
        while order2:
            Rw = E.random_element() # long Weierstrass coordinates
            order2 = Rw[1] == 0
        R = (Rw[0]*B, Rw[1]*B, Fq1) # Montgomery coordinates
        x3 = double_affine_x_only(R[0], A) # Montgomery coordinates
        Qw = 2*Rw
        Q = (Qw[0]*B, Qw[1]*B, Qw[2]) # Montgomery coordinates
        ok = x3 == Q[0]
        if not ok:
            print("Error double_affine_x_only")
            print("input: R[0]={}".format(R))
            print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
            print("got out: x3={}".format(x3))
        i = i+1
    print("test double_affine_x_only with B={}: {}".format(B, ok))
    return ok

def test_add_proj_xz_only(E, A, B):
    """
    test projective addition of points in Montgomery form and Montgomery X,Z coordinates
    E is in SageMath-compatible Montgomery form: y^2 = x^3 + A/B*x^2 + 1/B^2*x
    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    list_points = get_torsion_points(E, A, B)
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_l = [Fq(l) for l in [1, -1, 2, 3, 11]] # a set of scalars to test the projective coordinates (x,y)->(x*l,y*l,l)
    ok = True
    i0 = 0
    while ok and i0 < len(list_points):
        R = list_points[i0]
        i1 = 0
        while ok and i1 < len(list_points):
            S = list_points[i1]
            # the add_proj_xz_only function assumes that Z1*Z2 != 0 and X1*Z2 - X2*Z1 != 0
            # if R[0]/R[2] == S[0]/S[2], X4/Z4 is not defined with this formula.
            if R[0] == S[0] and R[2] != 0 and S[2] != 0:
                i1 = i1+1
                continue
            Rw = E(R[0]/B, R[1]/B, R[2]) # long Weierstrass coordinates
            Sw = E(S[0]/B, S[1]/B, S[2]) # long Weierstrass coordinates
            Tw = Rw - Sw # long Weierstrass coordinates
            T = (Tw[0] * B, Tw[1]*B, Tw[2])
            Qw = Rw + Sw
            Q = (Qw[0]*B, Qw[1]*B, Qw[2])
            j0 = 0
            while ok and j0 < len(list_l):
                lq = list_l[j0]
                Rp = (R[0]*lq, R[1]*lq, R[2]*lq) # projective coordinates
                j1 = 0
                while ok and j1 < len(list_l):
                    mq = list_l[j1]
                    Sp = (S[0]*mq, S[1]*mq, S[2]*mq) # projective coordinates
                    j2 = 0
                    while ok and j2 < len(list_l):
                        nq = list_l[j2]
                        Tp = (T[0]*nq, T[1]*nq, T[2]*nq) # projective coordinates
                        X3, Z3 = add_proj_xz_only(Rp[0], Rp[2], Sp[0], Sp[2], Tp[0], Tp[2])
                        ok = Q[0]*Z3 == X3
                        if not ok:
                            print("# Error add_proj_xz_only, inputs:")
                            print("  R={}".format(R))
                            print("  S={}".format(S))
                            print("  T={}".format(T))
                            print("  Rp=({}, {})\n  Sp=({}, {})\n  Tp=({}, {})".format(Rp[0], Rp[2], Sp[0], Sp[2], Tp[0], Tp[2]))
                            print("# expected output:\n x3={} #(and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
                            print("# got output:\n X3,Z3={}, {}".format(X3, Z3))
                            if Z3 != 0:
                                print("#  X3/Z3={}".format(X3/Z3))
                        j2 = j2+1
                    j1 = j1+1
                j0 = j0+1
            i1 = i1+1
        i0 = i0+1
    i = 0
    while ok and i < 10:
        distinct=False
        while not distinct:
            Rw = E.random_element() # long Weierstrass coordinates
            Sw = E.random_element() # long Weierstrass coordinates
            distinct = Rw[0] != Sw[0]
        Tw = Rw - Sw
        Qw = Rw + Sw # long Weierstrass coordinates
        Q = (Qw[0]*B, Qw[1]*B, Qw[2])
        R = (Rw[0]*B, Rw[1]*B, Fq1) # Montgomery coordinates
        S = (Sw[0]*B, Sw[1]*B, Fq1) # Montgomery coordinates
        T = (Tw[0]*B, Tw[1]*B, Fq1) # Montgomery coordinates
        j0 = 0
        while ok and j0 < len(list_l):
            lq = list_l[j0]
            Rp = (R[0]*lq, R[1]*lq, R[2]*lq) # projective coordinates
            j1 = 0
            while ok and j1 < len(list_l):
                mq = list_l[j1]
                Sp = (S[0]*mq, S[1]*mq, S[2]*mq) # projective coordinates
                j2 = 0
                while ok and j2 < len(list_l):
                    nq = list_l[j2]
                    Tp = (T[0]*nq, T[1]*nq, T[2]*nq) # projective coordinates
                    X3, Z3 = add_proj_xz_only(R[0], R[2], S[0], S[2], T[0], T[2]) # Montgomery coordinates
                    ok = X3 == Q[0]*Z3
                    if not ok:
                        print("# Error add_proj_xz_only, inputs:")
                        print("  R={}".format(R))
                        print("  S={}".format(S))
                        print("  T={}".format(T))
                        print("  Rp=({}, {})\n  Sp=({}, {})\n  Tp=({}, {})".format(Rp[0], Rp[2], Sp[0], Sp[2], Tp[0], Tp[2]))
                        print("# expected output:\n x3={} #(and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
                        print("# got output:\n X3,Z3={}, {}".format(X3, Z3))
                        if Z3 != 0:
                            print("#  X3/Z3={}".format(X3/Z3))
                    j2 = j2+1
                j1 = j1+1
            j0 = j0+1
        i = i+1
    print("test add_proj_xz_only with B={}: {}".format(B, ok))
    return ok

def test_mixed_add_proj_xz_only(E, A, B):
    """
    test mixed projective addition of points in Montgomery form and Montgomery X,Z coordinates
    E is in SageMath-compatible Montgomery form: y^2 = x^3 + A/B*x^2 + 1/B^2*x
    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    list_points = get_torsion_points(E, A, B)
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_l = [Fq(l) for l in [1, -1, 2, 3, 11]]
    ok = True
    i0 = 0
    while ok and i0 < len(list_points):
        R = list_points[i0]
        i1 = 0
        while ok and i1 < len(list_points):
            S = list_points[i1]
            # the mixed_add_proj_xz_only function assumes that X1*Z2 - X2*Z1 != 0, if Z1*Z2 != 0
            # if R[0]/R[2] == S[0]/S[2], X4/Z4 is not defined with this formula.
            if R[0] == S[0] and R[2] != 0 and S[2] != 0:
                i1 = i1+1
                continue
            Rw = E(R[0]/B, R[1]/B) # long Weierstrass coordinates
            Sw = E(S[0]/B, S[1]/B) # long Weierstrass coordinates
            Tw = Rw - Sw # long Weierstrass coordinates
            T = (Tw[0]*B, Tw[1]*B, Tw[2]*Fq1)
            Qw = Rw + Sw
            Q = (Qw[0]*B, Qw[1]*B, Qw[2])
            j0 = 0
            while ok and j0 < len(list_l):
                lq = list_l[j0]
                Rp = (R[0]*lq, R[1]*lq, R[2]*lq) # projective coordinates
                j1 = 0
                while ok and j1 < len(list_l):
                    mq = list_l[j1]
                    Sp = (S[0]*mq, S[1]*mq, S[2]*mq) # projective coordinates
                    X3, Z3 = mixed_add_proj_xz_only(Rp[0], Rp[2], Sp[0], Sp[2], T[0])
                    ok = Q[0]*Z3 == X3
                    if not ok:
                        print("Error mixed_add_proj_xz_only")
                        print("Inputs R={} S={} T={}".format(R, S, T))
                        print("Inputs Rp=({}, {}) Sp=({}, {}) T={}".format(Rp[0], Rp[2], Sp[0], Sp[2], T[0]))
                        print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
                        print("got out: X3,Z3={}, {}".format(X3, Z3))
                        if Z3 != 0:
                            print("got out: X3/Z3={}".format(X3/Z3))
                    j1 = j1+1
                j0 = j0+1
            i1 = i1+1
        i0 = i0+1
    i = 0
    while ok and i < 10:
        distinct=False
        while not distinct:
            Rw = E.random_element() # long Weierstrass coordinates
            Sw = E.random_element() # long Weierstrass coordinates
            distinct = Rw[0] != Sw[0]
        Tw = Rw - Sw
        Qw = Rw + Sw # long Weierstrass coordinates
        Q = (Qw[0]*B, Qw[1]*B, Qw[2])
        R = (Rw[0]*B, Rw[1]*B, Fq1) # Montgomery coordinates
        S = (Sw[0]*B, Sw[1]*B, Fq1) # Montgomery coordinates
        T = (Tw[0]*B, Tw[1]*B, Tw[2]) # Montgomery coordinates
        j0 = 0
        while ok and j0 < len(list_l):
            lq = list_l[j0]
            Rp = (R[0]*lq, R[1]*lq, R[2]*lq) # projective coordinates
            j1 = 0
            while ok and j1 < len(list_l):
                mq = list_l[j1]
                Sp = (S[0]*mq, S[1]*mq, S[2]*mq) # projective coordinates
                X3, Z3 = mixed_add_proj_xz_only(R[0], R[2], S[0], S[2], T[0]) # Montgomery coordinates
                ok = X3 == Q[0]*Z3
                if not ok:
                    print("Error mixed_add_proj_xz_only")
                    print("Inputs R={} S={} T={}".format(R, S, T))
                    print("Inputs Rp=({}, {}) Sp=({}, {}) T={}".format(Rp[0], Rp[2], Sp[0], Sp[2], T[0]))
                    print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
                    print("got out: X3,Z3={}, {}".format(X3, Z3))
                    if Z3 != 0:
                        print("got out: X3/Z3={}".format(X3/Z3))
                j1 = j1+1
            j0 = j0+1
        i = i+1
    print("test mixed_add_proj_xz_only with B={}: {}".format(B, ok))
    return ok

def test_double_proj_xz_only(E, A, B):
    """
    test doubling of points in projective Montgomery form and Montgomery (X,Z)-only coordinates
    E is in SageMath-compatible Montgomery form: y^2 = x^3 + A/B*x^2 + 1/B^2*x
    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    list_points = get_torsion_points(E, A, B)
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_l = [Fq(l) for l in [1, -1, 2, 3, 11]]
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        # the double_proj_xz_only function does not assume anything on the data
        # if the point is a 2-torsion point (0,0,1) or (a,0,1) it will return (0,0) in (X,Z)
        Rw = E(R[0]/B, R[1]/B, R[2]) # long Weierstrass coordinates
        Qw = 2*Rw # long Weierstrass coordinates
        Q = (Qw[0]*B, Qw[1]*B, Qw[2])
        X3, Z3 = double_proj_xz_only(R[0], R[2], A)
        ok = Q[0]*Z3 == X3
        if not ok:
            print("Error double_proj_xz_only")
            print("input: R[0]={} R[2]={}".format(R[0], R[2]))
            print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
            print("got out: X3,Z3={}, {}".format(X3, Z3))
            if Z3 != 0:
                print("got out: X3/Z3={}".format(X3/Z3))
        i = i+1
    i = 0
    while ok and i < 10:
        Rw = E.random_element() # long affine Weierstrass coordinates (X, Y, 1)
        R = (Rw[0]*B, Rw[1]*B, Fq1) # affine Montgomery coordinates
        Qw = 2*Rw
        Q = (Qw[0]*B, Qw[1]*B, Fq1) # affine Montgomery coordinates
        j = 0
        while ok and j < len(list_l):
            l = list_l[j]
            X3, Z3 = double_proj_xz_only(R[0]*l, R[2]*l, A) # Montgomery coordinates
            ok = X3 == Q[0]*Z3
            if not ok:
                print("Error double_proj_xz_only i={} j={}".format(i, j))
                print("input: R[0]={} R[2]={}".format(R[0]*l, R[2]*l))
                print("exp out: x3={} (and y3={}, z3={})".format(Q[0], Q[1], Q[2]))
                print("got out: X3,Z3={}, {}".format(X3, Z3))
                if Z3 != 0:
                    print("got out: X3/Z3={}".format(X3/Z3))
            j = j+1
        i = i+1
    print("test double_proj_xz_only with B={}: {}".format(B, ok))
    return ok

def test_montgomery_ladder_affine_x_only(E, A, B, prime_r, cofactor):
    """
    Test the Montgomery ladder in affine x-only coordinates

    montgomery_ladder_affine_x_only(m, P, A)
    INPUT:
    - `E`: elliptic curve in SageMath compatible format y^2 = x^3 + A/B*x^2 + 1/B^2*x
    - `A`: curve coefficient
    - `B`: curve coefficient
    - `prime_r`: prime order of the subgroup
    - `cofactor`: cofactor of the subgroup, so that E.order() = p+1-t = cofactor * prime_r

    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    """
    # check that E is in SageMath compatible Montgomery form
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    i = 0
    ok = True
    order = E.order()
    while ok and i < 10:
        Pw = E.random_element()
        P = (Pw[0] * B, Pw[1] * B, Pw[2])
        m = randrange(order)
        mPw = m*Pw
        mP = (mPw[0] * B, mPw[1] * B, mPw[2])
        x3 = montgomery_ladder_affine_x_only(m, P, A)
        ok = mP[0] == x3
        if not ok:
            print("error montgomery_ladder_affine_x_only")
            print("input m = {} P = {} A = {}".format(m, P, A))
            print("expected output m*P = {}".format(mP))
            print("got output x3 = {}".format(x3))
        i = i+1
    print("test montgomery_ladder_affine_x_only with B={}: {}".format(B, ok))
    return ok

def test_montgomery_ladder_proj_xz_only(E, A, B, prime_r, cofactor):
    """
    Test the Montgomery ladder in projective X,Z-only coordinates

    montgomery_ladder_proj_xz_only(m, P, A)
    INPUT:
    - `E`: elliptic curve in SageMath compatible format y^2 = x^3 + A/B*x^2 + 1/B^2*x
    - `A`: curve coefficient
    - `B`: curve coefficient
    - `prime_r`: prime order of the subgroup
    - `cofactor`: cofactor of the subgroup, so that E.order() = p+1-t = cofactor * prime_r

    map from E to EM : B*y^2 = x^3 + A*x^2 + x is (x,y) -> (x*B, y*B, 1)
    """
    # check that E is in SageMath compatible Montgomery form
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    i = 0
    ok = True
    order = E.order()
    m = 0
    while ok and m < 10:
        Pw = E.random_element()
        P = (Pw[0] * B, Pw[1] * B, 1) # affine coordinates
        mPw = m*Pw
        mP = (mPw[0] * B, mPw[1] * B, mPw[2])
        X3, Z3 = montgomery_ladder_proj_xz_only(m, P, A)
        ok = mP[0] * Z3 == X3
        if not ok:
            print("error montgomery_ladder_proj_xz_only")
            print("input m = {} P = {} A = {}".format(m, P, A))
            print("expected output m*P = {}".format(mP))
            print("got output X3, Z3 = {}, {}".format(X3, Z3))
            if Z3 != 0:
                print("got X3/Z3 = {}".format(X3/Z3))
        m = m+1
    i = 0
    while ok and i < 10:
        Pw = E.random_element()
        P = (Pw[0] * B, Pw[1] * B, Pw[2])
        m = randrange(order)
        mPw = m*Pw
        mP = (mPw[0] * B, mPw[1] * B, mPw[2])
        X3, Z3 = montgomery_ladder_proj_xz_only(m, P, A)
        ok = mP[0] * Z3 == X3
        if not ok:
            print("error montgomery_ladder_proj_xz_only")
            print("input m = {} P = {} A = {}".format(m, P, A))
            print("expected output m*P = {}".format(mP))
            print("got output X3, Z3 = {}, {}".format(X3, Z3))
            if Z3 != 0:
                print("got X3/Z3 = {}".format(X3/Z3))
        i = i+1
    print("test montgomery_ladder_proj_xz_only with B={}: {}".format(B, ok))
    return ok

################################################################################
## TESTS GIVEN AS EXAMPLES
################################################################################
def get_torsion_points(E, A, B):
    """
    Return the F-rational 2-torsion and 4-torsion points of E in Montgomery form
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    Fp = A.parent()
    Fp0 = Fp(0)
    Fp1 = Fp(1)
    # find the 2-torsion points
    # whenever B = 1 or not, the formulas are the same
    P = (Fp0, Fp0, Fp1)
    list_points = [P]
    if (A**2-4).is_square():
        alpha = (A**2-4).sqrt()
        P2 = ((-A + alpha)/2, Fp0, Fp1)
        P2_ = ((-A - alpha)/2, Fp0, Fp1)
        list_points.append(P2)
        list_points.append(P2_)
    # find the 4-torsion points
    if ((A+2)*B).is_square():
        delta = ((A+2)/B).sqrt()
        Q4 = (1, delta, Fp1)
        Q4_ = (1, -delta, Fp1)
        list_points.append(Q4)
        list_points.append(Q4_)
    if ((A-2)*B).is_square():
        gamma = ((A-2)/B).sqrt()
        R4 = (-1, gamma, Fp1)
        R4_ = (-1, -gamma, Fp1)
        list_points.append(R4)
        list_points.append(R4_)
    return list_points

def test_add_affine_montgomery(E, A, B):
    """
    test affine addition of points in Montgomery form
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    list_points = get_torsion_points(E, A, B)
    Fq = E.base_field()
    Fq1 = Fq(1)
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        for S in list_points:
            res = add_affine_montgomery(R, S, A, B)
            Rw = E(R[0]/B, R[1]/B) # long Weierstrass coordinates
            Sw = E(S[0]/B, S[1]/B) # long Weierstrass coordinates
            exp = Rw+Sw # long Weierstrass coordinates
            # special case for the point at infinity
            if (exp == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (exp[0]*B, exp[1]*B)
                ok = ok and res[0] == exp[0] and res[1] == exp[1] and res[2] == 1
            if not ok:
                print("input: R={}".format(R))
                print("input: S={}".format(S))
                print("exp out: {}".format(exp))
                print("got out: {}".format(res))
        i = i+1
    i = 0
    while ok and i < 10:
        Rw = E.random_element() # long Weierstrass coordinates
        Sw = E.random_element() # long Weierstrass coordinates
        R = (Rw[0]*B, Rw[1]*B, Fq1) # Montgomery coordinates
        S = (Sw[0]*B, Sw[1]*B, Fq1) # Montgomery coordinates
        res = add_affine_montgomery(R, S, A, B) # Montgomery coordinates
        exp = Rw+Sw # long Weierstrass coordinates
        if (exp == E(0)):
            ok = ok and res[2] == 0
        else:
            exp = (exp[0]*B, exp[1]*B)
            ok = ok and res[0] == exp[0] and res[1] == exp[1] and res[2] == 1
        if not ok:
            print("input: R={}".format(R))
            print("input: S={}".format(S))
            print("exp out: {}".format(exp))
            print("got out: {}".format(res))
        i = i+1
    print("test add_affine_montgomery with B={}: {}".format(B, ok))
    return ok

def test_dbl_affine_montgomery(E, A, B):
    """
    test affine doubling of points in Montgomery form
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_points = get_torsion_points(E, A, B)
    ok = True
    i = 0
    while i < len(list_points):
        R = list_points[i]
        res = dbl_affine_montgomery(R, A, B)
        Rw = E(R[0]/B, R[1]/B)
        exp = 2*Rw
        # special case for the point at infinity
        if (exp == E(0)):
            ok = ok and res[2] == 0
        else:
            exp = (exp[0]*B, exp[1]*B)
            ok = ok and res[0] == exp[0] and res[1] == exp[1] and res[2] == 1
        if not ok:
            print("input: R={}".format(R))
            print("exp out: {}".format(exp))
            print("got out: {}".format(res))
        i = i+1
    i = 0
    while ok and i < 10:
        Rw = E.random_element()
        while Rw[2] == 0:
            Rw = E.random_element()
        R = (Rw[0]*B, Rw[1]*B, Fq1)
        res = dbl_affine_montgomery(R, A, B)
        exp = 2*Rw
        # special case for the point at infinity
        if (exp == E(0)):
            ok = ok and res[2] == 0
        else:
            exp = (exp[0]*B, exp[1]*B)
            ok = ok and res[0] == exp[0] and res[1] == exp[1] and res[2] == 1
        if not ok:
            print("test {} wrong".format(i))
            print("input: R={}".format(R))
            print("exp out: {}".format(exp))
            print("got out: {}".format(res))
        i = i+1
    print("test dbl_affine_montgomery with B={}: {}".format(B, ok))
    return ok

def test_mixed_add_projective_montgomery(E, A, B):
    """
    test mixed addition of points (affine and projective) on E in Montgomery form
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_points = get_torsion_points(E, A, B)
    list_l = [1, -1, 2, 3, 11]
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        Rw = E(R[0]/B, R[1]/B) # long Weierstrass coordinates
        j = 0
        while ok and j < len(list_l):
            l = list_l[j]
            lq = Fq(l)
            Rp = (R[0]*lq, R[1]*lq, lq) # projective coordinates
            for S in list_points:
                Sw = E(S[0]/B, S[1]/B) # long Weierstrass coordinates
                exp = Rw+Sw # long Weierstrass coordinates
                res = mixed_add_proj_montgomery(Rp, S, A, B)
                # special case for the point at infinity
                if (exp == E(0)):
                    ok = ok and res[2] == 0
                else:
                    exp = (exp[0]*B, exp[1]*B)
                    ok = ok and res[0] == exp[0]*res[2] and res[1] == exp[1]*res[2]
                if not ok:
                    print("input: R={}".format(R))
                    print("input: Rp={}".format(Rp))
                    print("input: S={}".format(S))
                    print("exp out: {}".format(exp))
                    print("got out: {}".format(res))
                    print("got out: {}".format((res[0]/res[2], res[1]/res[2])))
                    X3, Y3, Z3 = res
                    point_on_curve = Y3**2*Z3 * B == X3**3 + A * X3**2 * Z3 + X3 * Z3**2
                    print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    i = 0
    while ok and i < 10:
        Rw = E.random_element() # long Weierstrass coordinates
        Sw = E.random_element() # long Weierstrass coordinates
        R = (Rw[0]*B, Rw[1]*B, Fq1) # Montgomery coordinates
        S = (Sw[0]*B, Sw[1]*B, Fq1) # Montgomery coordinates
        expw = Rw+Sw # long Weierstrass coordinates
        j = 0
        while ok and j < len(list_l):
            l = list_l[j]
            Rp = (R[0]*l, R[1]*l, l)
            res = mixed_add_proj_montgomery(Rp, S, A, B) # Montgomery coordinates
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0]*B, expw[1]*B)
                ok = ok and res[0] == exp[0]*res[2] and res[1] == exp[1]*res[2]
            if not ok:
                print("input: R={}".format(R))
                print("input: Rp={}".format(Rp))
                print("input: S={}".format(S))
                print("exp out: {}".format(exp))
                print("got out: {}".format(res))
                print("got out: {}".format((res[0]/res[2], res[1]/res[2])))
                X3, Y3, Z3 = res
                point_on_curve = Y3**2*Z3 * B == X3**3 + A * X3**2 * Z3 + X3 * Z3**2
                print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    print("test mixed_add_proj_montgomery with B={}: {}".format(B, ok))
    return ok

def test_dbl_projective_montgomery(E, A, B):
    """
    test doubling of a point in projective coordinates on E in Montgomery form
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    list_points = get_torsion_points(E, A, B)
    list_l = [1, -1, 2, 3, 11]
    Fp = A.parent()
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        Rw = E(R[0]/B, R[1]/B) # long Weierstrass coordinates
        expw = 2*Rw # long Weierstrass coordinates
        j = 0
        while ok and j < len(list_l):
            l = list_l[j]
            lq = Fp(l)
            Rp = (R[0]*lq, R[1]*lq, lq) # projective coordinates
            res = dbl_proj_montgomery(Rp, A, B)
            # special case for the point at infinity
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0]*B, expw[1]*B)
                ok = ok and res[0] == exp[0]*res[2] and res[1] == exp[1]*res[2]
            if not ok:
                print("input: R={}".format(R))
                print("input: Rp={}".format(Rp))
                print("exp out: {}".format(exp))
                print("got out: {}".format(res))
                print("got out: {}".format((res[0]/res[2], res[1]/res[2])))
                X3, Y3, Z3 = res
                point_on_curve = Y3**2*Z3 * B == X3**3 + A * X3**2 * Z3 + X3 * Z3**2
                print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    i = 0
    while ok and i < 10:
        Rw = E.random_element() # long Weierstrass coordinates
        R = (Rw[0]*B, Rw[1]*B) # Montgomery coordinates
        expw = 2*Rw # long Weierstrass coordinates
        j = 0
        while ok and j < len(list_l):
            l = list_l[j]
            Rp = (R[0]*l, R[1]*l, l)
            res = dbl_proj_montgomery(Rp, A, B) # Montgomery coordinates
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0]*B, expw[1]*B)
                ok = ok and res[0] == exp[0]*res[2] and res[1] == exp[1]*res[2]
            if not ok:
                print("input: R={}".format(R))
                print("input: Rp={}".format(Rp))
                print("exp out: {}".format(exp))
                print("got out: {}".format(res))
                print("got out: {}".format((res[0]/res[2], res[1]/res[2])))
                X3, Y3, Z3 = res
                point_on_curve = Y3**2*Z3 * B == X3**3 + A * X3**2 * Z3 + X3 * Z3**2
                print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    print("test dbl_proj_montgomery with B={}: {}".format(B, ok))
    return ok

def test_double_and_add_proj_montgomery(E, A, B):
    """
    test double_and_add_proj_montgomery(m, P, A, B)
    """
    assert E.a1() == 0 and E.a3() == 0 and E.a6() == 0 and B**2*E.a4() == 1 and B*E.a2() == A
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_l = [1, -1, 2, -2, 3, -3, 7, -7, 8, -8, 11, -11]
    i = 0
    ok = True
    while ok and i < 10:
        Pw = E.random_element()
        P = (Pw[0]*B, Pw[1]*B, Fq1)
        j = 0
        while ok and j < len(list_l):
            m = list_l[j]
            expw = m * Pw
            res = double_and_add_proj_montgomery(m, P, A, B)
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0]*B, expw[1]*B)
                ok = ok and res[0] == exp[0]*res[2] and res[1] == exp[1]*res[2]
                if not ok:
                    print("input: P={}".format(P))
                    print("input: m={}".format(m))
                    print("exp out: {}".format(exp))
                    print("got out: {}".format(res))
                    print("got out: {}".format((res[0]/res[2], res[1]/res[2])))
                    X3, Y3, Z3 = res
                    point_on_curve = Y3**2*Z3 * B == X3**3 + A * X3**2 * Z3 + X3 * Z3**2
                    print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    print("test double_and_add_proj_montgomery with B={}: {}".format(B, ok))
    return ok


if __name__ == "__main__":
    p = ZZ(2**255-19)
    Fp = GF(p)
    # y^2 = x^3 + 486662*x^2 + x is one form of this curve
    # Montgomery form is
    A = Fp(486662)
    B = Fp(1)
    EM = EllipticCurve([0, A/B, 0, 1/B**2, 0]) # it will throw an exception if EM is singular

    # two options for computing the order: method EM.order() or with the trace
    orderE = EM.order()
    assert orderE % 4 == 0
    r = orderE // 8 # actually the curve order is multiple of 8
    cofactor = ZZ(8)
    assert r.is_prime() # the curve order is 8 times a large prime of 253 bits
    tr = EM.trace_of_frobenius() # alternative
    orderEtr = p + 1 - tr
    assert orderE == orderEtr
    # base point:
    x0 = Fp(9)
    y02 = x0**3 + A*x0**2 + B*x0
    y0 = y02.sqrt()
    P = EM((x0, y0))
    assert r*P == EM(0)
    assert not EM.is_supersingular() # the curve is not supersingular
    assert gcd(tr, p) == 1 # alternative: the trace is coprime to the characteristic p

    test_add_affine_montgomery(EM, A, B)
    test_dbl_affine_montgomery(EM, A, B)
    test_mixed_add_projective_montgomery(EM, A, B)
    test_dbl_projective_montgomery(EM, A, B)
    test_double_and_add_proj_montgomery(EM, A, B)

    test_add_affine_x_only(EM, A, B)
    test_double_affine_x_only(EM, A, B)
    test_add_proj_xz_only(EM, A, B)
    test_mixed_add_proj_xz_only(EM, A, B)
    test_double_proj_xz_only(EM, A, B)
    test_montgomery_ladder_affine_x_only(EM, A, B, r, cofactor)
    test_montgomery_ladder_proj_xz_only(EM, A, B, r, cofactor)
    B = Fp(2)
    while B.is_square():
        B = B+1
    print("B= {} is not a square in Fp".format(B))
    EB = EllipticCurve([0, A/B, 0, 1/B**2, 0])
    assert EB.is_isomorphic(EM.quadratic_twist())
    #EMt = EM.quadratic_twist()
    #phi = EB.isomorphism_to(EMt)
    #print(phi)
    test_add_affine_montgomery(EB, A, B)
    test_dbl_affine_montgomery(EB, A, B)
    test_mixed_add_projective_montgomery(EB, A, B)    
    test_dbl_projective_montgomery(EB, A, B)
    test_double_and_add_proj_montgomery(EB, A, B)

    test_add_affine_x_only(EB, A, B)
    test_double_affine_x_only(EB, A, B)
    test_add_proj_xz_only(EB, A, B)
    test_mixed_add_proj_xz_only(EB, A, B)
    test_double_proj_xz_only(EB, A, B)
    test_montgomery_ladder_affine_x_only(EB, A, B, r, cofactor)
    test_montgomery_ladder_proj_xz_only(EB, A, B, r, cofactor)