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_char2(R, a_2, a_6):
    """
    Doubling on a binary curve E2: y^2 + xy = x^3 + a_2*x^2 + a_6

    INPUT:
    - `R`: point in affine coordinates (x, y, 1) or (0, y, 0) (point at infinity)
    - `a_2`: curve coefficient on E2
    - `a_6`: curve coefficient on E2

    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
    https://www.hyperelliptic.org/EFD/g12o/auto-shortw-affine.html
    """
    Fq = a_2.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_char2 the input point R is not in affine coordinates (x1, y1, 1): R = {}".format(R))
    if x1 == 0:
        return (Fq0, Fq1, Fq0) # point at infinity
    l = x1 + y1/x1
    x3 = l**2 + l + a_2
    y3 = l * (x1 + x3) + x3 + y1
    return (x3, y3, Fq1)

def add_affine_char2(P, R, a_2, a_6):
    """
    Doubling on a binary curve E2: y^2 + xy = x^3 + a_2*x^2 + a_6

    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_2`: curve coefficient on E2
    - `a_6`: curve coefficient on E2

    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
    https://www.hyperelliptic.org/EFD/g12o/auto-shortw-affine.html
    """
    Fq = a_2.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_char2 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_char2 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 == x2+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_char2(P, a_2, a_6)
    # here we can assume that (y1 != 0 or y2 != 0) or (x1 != x2)
    l = (y2 + y1)/(x1 + x2)
    x3 = l**2 + l + x1 + x2 + a_2
    y3 = l * (x1 + x3) + x3 + y1
    return (x3, y3, 1)

def dbl_proj_char2(R, a_2, a_6):
    """
    Doubling on on a binary curve E: Y^2*Z + XYZ = X^3 + a_2*X^2Z + a_6*Z^3

    INPUT:
    - `R`: point in projective coordinates (X,Y,Z), if Z = 0 the point is at infinity
    - `a_2`: curve coefficient on E2
    - `a_6`: curve coefficient on E2

    RETURN: Q(X3, Y3, Z3) = 2*R in projective coordinates
    """
    Fq = a_2.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)
    A = X1**2
    B = A+Y1*Z1
    C = X1*Z1
    BC = B+C
    D = C**2
    E = B*BC+a_2*D
    X3 = C*E
    Y3 = BC*E+A**2*C
    Z3 = C*D
    return (X3, Y3, Z3)

def mixed_add_proj_char2(R, P, a_2, a_6):
    """
    Mixed Addition on a binary curve E: Y^2*Z + XYZ = X^3 + a_2*X^2Z + a_6*Z^3

    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_2`: curve coefficient on E2
    - `a_6`: curve coefficient on E2

    RETURN: Q(X3, Y3, Z3) = R+P in projective coordinates
    """
    Fq = a_2.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_char2 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 == y2*Z1 + X1): # x2 == X1/Z1 and Y1/Z1 == x2+y2 <=> Y1 == (x2+y2)*Z1 <=> y2*Z1 + X1, 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_char2(P, a_2, a_6)
    A = Y1+Z1*y2
    B = X1+Z1*x2
    AB = A+B
    C = B**2
    E = B*C
    F = (A*AB+a_2*C)*Z1+E
    X3 = B*F
    Y3 = C*(A*X1+B*Y1)+AB*F
    Z3 = E*Z1
    return (X3, Y3, Z3)

def double_and_add_proj_char2(m, P, a_2, a_6):
    """
    scalar multiplication [m]*P in binary curve E: Y^2*Z + XYZ = X^3 + a_2*X^2Z + a_6*Z^3

    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_2.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[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_char2(R, a_2, a_6)
        if b[i] == 1:
            R = mixed_add_proj_char2(R, P, a_2, a_6)
    return R

def add_affine_x_only_char2(x1, x2, x4):
    """
    x-only addition for the Montgomery ladder on E: y^2 + xy = x^3 + a_2*x^2 + a_6

    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)
    """
    x3 = x4+(x1*x2)/(x1+x2)**2
    return x3

def double_affine_x_only_char2(x1, a_6):
    """
    x-only doubling for the Montgomery ladder on E: y^2 + xy = x^3 + a_2*x^2 + a_6

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

    RETURN: x3 the x-coordinate of [2]P
    """
    # TODO your code here, and change the return statement
    x3 = x1**2 + a_6/x1**2
    return x3

def add_proj_xz_only_char2(X1, Z1, X2, Z2, X4, Z4):
    """
    Addition in Projective X,Z-only on a binary curve E: Y^2*Z + XYZ = X^3 + a_2*X^2Z + a_6*Z^3

    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)
    A = X1*X2
    B = Z1*Z2
    X3 = Z4*(A**2+a_6*B**2)
    Z3 = X4*((X1+Z1)*(X2+Z2)-A-B)**2
    return X3, Z3

def mixed_add_proj_xz_only_char2(X1, Z1, X2, Z2, x4):
    """
    Mixed Addition in Projective X,Z-only on a binary curve E: Y^2*Z + XYZ = X^3 + a_2*X^2Z + a_6*Z^3

    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
    https://www.hyperelliptic.org/EFD/g12o/auto-shortw-xz.html#diffadd-mdadd-2003-s

    Exceptional cases:
    if X1/Z1 = X2/Z2 and Z1*Z2 != 0 <=> X1*Z2 - X2*Z1 = 0 then it should raise an exception
    """
    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)
    A = X1*Z2
    B = X2*Z1
    Z3 = (A+B)**2
    X3 = x4*Z3+A*B
    return X3, Z3

def double_proj_xz_only_char2(X1, Z1, a_6):
    """
    Doubling in Projective X,Z-only on a binary curve E: Y^2*Z + XYZ = X^3 + a_2*X^2Z + a_6*Z^3

    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

    if the point is a 2-torsion point it will return (0,0) in (X,Z)
    """
    if Z1 == 0 or X1 == 0:
        return (0, 0)
    XX1 = X1**2
    ZZ1 = Z1**2
    X3 = (XX1+sqrt(a_6)*ZZ1)**2
    Z3 = XX1*ZZ1
    return X3, Z3

# Montgomery ladder in affine and projective coordinates
def montgomery_ladder_affine_x_only_char2(m, P, a_6):
    """
    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_char2(xP, a_6)
    for i in range(len(mbits)-2, -1, -1):
        bi = mbits[i]
        if bi == 0:
            (x0, x1) = double_affine_x_only_char2(x0, a_6), add_affine_x_only_char2(x0, x1, xP)
        else:
            (x0, x1) = add_affine_x_only_char2(x0, x1, xP), double_affine_x_only_char2(x1, a_6)
    return x0

def montgomery_ladder_proj_xz_only_char2(m, P, a_2, a_6, E=None):
    """
    Montgomery ladder in projective 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]
    ZP = P[2]
    assert ZP == 1
    X0 = (XP, ZP)
    X1 = double_proj_xz_only_char2(XP, ZP, a_6)
    for i in range(len(mbits)-2, -1, -1):
        bi = mbits[i]
        if bi == 0:
            (X0, X1) = double_proj_xz_only_char2(X0[0], X0[1], a_6), mixed_add_proj_xz_only_char2(X0[0], X0[1], X1[0], X1[1], XP)
        else:
            (X0, X1) = mixed_add_proj_xz_only_char2(X0[0], X0[1], X1[0], X1[1], XP), double_proj_xz_only_char2(X1[0], X1[1], a_6)
    return X0

################################################################################
## TESTS FOR THE HAND-IN
################################################################################
def test_add_affine_x_only_char2(E, a_2, a_6):
    """
    test affine addition of points in x-only coordinates
    It is assumed that the inputs are consistent, that is the points are distincts, non-opposite, and non-zero
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    # first, test addition of a torsion point with another distinct random point
    list_points = get_torsion_points(E, a_2, a_6)
    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], R[1], R[2])  # long Weierstrass coordinates
        S = (Sw[0], Sw[1], Sw[2]) # Montgomery coordinates
        Tw = Rw - Sw                  # long Weierstrass coordinates
        T = (Tw[0], Tw[1], Tw[2])
        x4 = T[0]
        x3 = add_affine_x_only_char2(R[0], S[0], x4)
        Qw = Rw + Sw
        Q = (Qw[0], Qw[1], 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], Rw[1], Fq1) # Montgomery coordinates
        S = (Sw[0], Sw[1], Fq1) # Montgomery coordinates
        T = (Tw[0], Tw[1], Fq1) # Montgomery coordinates
        x3 = add_affine_x_only_char2(R[0], S[0], T[0]) # Montgomery coordinates
        Qw = Rw+Sw # long Weierstrass coordinates
        Q = (Qw[0], Qw[1], 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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_double_affine_x_only_char2(E, a_2, a_6):
    """
    test affine doubling of points in x-only coordinates in characteristic 2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    list_points = get_torsion_points(E, a_2, a_6)
    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_char2 function assumes that x1 is non-zero
        if R[0] == 0:
            i=i+1
            continue # it will skip the rest of the iteration and directly go to the next i
        Rw = E(R[0], R[1]) # long Weierstrass coordinates
        Qw = 2*Rw # long Weierstrass coordinates
        Q = (Qw[0], Qw[1], Qw[2])
        x3 = double_affine_x_only_char2(R[0], a_2)
        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], Rw[1], Fq1) # Montgomery coordinates
        x3 = double_affine_x_only_char2(R[0], a_2) # Montgomery coordinates
        Qw = 2*Rw
        Q = (Qw[0], Qw[1], Qw[2]) # Montgomery coordinates
        ok = x3 == Q[0]
        if not ok:
            print("Error double_affine_x_only_char2")
            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 a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_add_proj_xz_only_char2(E, a_2, a_6):
    """
    test projective addition of points in X,Z coordinates in characteristic 2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    list_points = get_torsion_points(E, a_2, a_6)
    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], R[1], R[2]) # long Weierstrass coordinates
            Sw = E(S[0], S[1], S[2]) # long Weierstrass coordinates
            Tw = Rw - Sw # long Weierstrass coordinates
            T = (Tw[0], Tw[1], Tw[2])
            Qw = Rw + Sw
            Q = (Qw[0], Qw[1], 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_char2, 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], Qw[1], Qw[2])
        R = (Rw[0], Rw[1], Fq1) # Montgomery coordinates
        S = (Sw[0], Sw[1], Fq1) # Montgomery coordinates
        T = (Tw[0], Tw[1], 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_char2(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_char2, 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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_mixed_add_proj_xz_only_char2(E, a_2, a_6):
    """
    test mixed projective addition of points in X,Z coordinates in characteristic 2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    list_points = get_torsion_points(E, a_2, a_6)
    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], R[1]) # long Weierstrass coordinates
            Sw = E(S[0], S[1]) # long Weierstrass coordinates
            Tw = Rw - Sw # long Weierstrass coordinates
            T = (Tw[0], Tw[1], Tw[2]*Fq1)
            Qw = Rw + Sw
            Q = (Qw[0], Qw[1], 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_char2(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_char2")
                        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], Qw[1], Qw[2])
        R = (Rw[0], Rw[1], Fq1) # Montgomery coordinates
        S = (Sw[0], Sw[1], Fq1) # Montgomery coordinates
        T = (Tw[0], Tw[1], 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_char2(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_char2")
                    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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_double_proj_xz_only_char2(E, a_2, a_6):
    """
    test doubling½ of points in X,Z coordinates in characteristic 2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    list_points = get_torsion_points(E, a_2, a_6)
    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], R[1], R[2]) # long Weierstrass coordinates
        Qw = 2*Rw # long Weierstrass coordinates
        Q = (Qw[0], Qw[1], Qw[2])
        X3, Z3 = double_proj_xz_only_char2(R[0], R[2], a_6)
        ok = Q[0]*Z3 == X3
        if not ok:
            print("Error double_proj_xz_only_char2")
            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 Weierstrass coordinates
        R = (Rw[0], Rw[1], Rw[2]) # Montgomery coordinates
        X3, Z3 = double_proj_xz_only_char2(R[0], R[2], a_6) # Montgomery coordinates
        Qw = 2*Rw
        Q = (Qw[0], Qw[1], Qw[2]) # Montgomery coordinates
        ok = X3 == Q[0]*Z3
        if not ok:
            print("Error double_proj_xz_only_char2")
            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
    print("test double_proj_xz_only_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_montgomery_ladder_affine_x_only_char2(E, a_2, a_6, prime_r, cofactor):
    """
    Test the Montgomery ladder in affine x-only coordinates

    montgomery_ladder_affine_x_only_char2(m, P, A)
    INPUT:
    - `E`: elliptic curve in SageMath compatible format y^2 + xy = x^3 + a_2*x^2 + a_6
    - `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
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    i = 0
    ok = True
    order = E.order()
    while ok and i < 10:
        Pw = E.random_element()
        P = (Pw[0], Pw[1], Pw[2])
        m = randrange(order)
        mPw = m*Pw
        mP = (mPw[0], mPw[1], mPw[2])
        x3 = montgomery_ladder_affine_x_only_char2(m, P, a_6)
        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 a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_montgomery_ladder_proj_xz_only_char2(E, a_2, a_6, prime_r, cofactor):
    """
    Test the Montgomery ladder in projective X,Z-only coordinates

    test_montgomery_ladder_proj_xz_only_char2(m, P, a_6)
    INPUT:
    - `E`: elliptic curve in SageMath compatible format y^2 + xy = x^3 + a_2*x^2 + a_6
    - `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
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    i = 0
    ok = True
    order = E.order()
    m = 0
    while ok and m < 10:
        Pw = cofactor * E.random_element()
        P = (Pw[0], Pw[1], Pw[2])
        mPw = m*Pw
        mP = (mPw[0], mPw[1], mPw[2])
        X3, Z3 = montgomery_ladder_proj_xz_only_char2(m, P, a_2, a_6, E)
        ok = mP[0] * Z3 == X3
        if not ok:
            print("error montgomery_ladder_proj_xz_only_char2")
            print("input m = {} P = {} a_6 = {}".format(m, P, a_6.integer_representation()))
            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 = cofactor * E.random_element()
        P = (Pw[0], Pw[1], Pw[2])
        m = randrange(order)
        mPw = m*Pw
        mP = (mPw[0], mPw[1], mPw[2])
        X3, Z3 = montgomery_ladder_proj_xz_only_char2(m, P, a_2, a_6, E)
        ok = mP[0] * Z3 == X3
        if not ok:
            print("error montgomery_ladder_proj_xz_only_char2")
            print("input m = {} P = {} a_6 = {}".format(m, P, a_6.integer_representation()))
            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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

################################################################################
## TESTS GIVEN AS EXAMPLES
################################################################################
def get_torsion_points(E, a_2, a_6):
    """
    Return the 2-torsion points
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    F = a_2.parent()
    Fq0 = F(0)
    Fq1 = F(1)
    # find the 2-torsion points
    P = (Fq0, sqrt(a_6), Fq1)
    list_points = [P]
    return list_points

def test_add_affine_char2(E, a_2, a_6):
    """
    test affine addition of points in characteristic 2
    """
    list_points = get_torsion_points(E, a_2, a_6)
    F = E.base_field()
    Fq1 = F(1)
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        for S in list_points:
            res = add_affine_char2(R, S, a_2, a_6)
            Rw = E(R[0], R[1])
            Sw = E(S[0], S[1])
            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], exp[1])
                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], Rw[1], Fq1) # Montgomery coordinates
        S = (Sw[0], Sw[1], Fq1) # Montgomery coordinates
        res = add_affine_char2(R, S, a_2, a_6) # Montgomery coordinates
        exp = Rw+Sw # long Weierstrass coordinates
        if (exp == E(0)):
            ok = ok and res[2] == 0
        else:
            exp = (exp[0], exp[1])
            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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_dbl_affine_char2(E, a_2, a_6):
    """
    test affine doubling of points in characteristic 2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    F = E.base_field()
    Fq1 = F(1)
    list_points = get_torsion_points(E, a_2, a_6)
    ok = True
    i = 0
    while i < len(list_points):
        R = list_points[i]
        res = dbl_affine_char2(R, a_2, a_6)
        Rw = E(R[0], R[1])
        exp = 2*Rw
        # special case for the point at infinity
        if (exp == E(0)):
            ok = ok and res[2] == 0
        else:
            exp = (exp[0], exp[1])
            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], Rw[1], Fq1)
        res = dbl_affine_char2(R, a_2, a_6)
        exp = 2*Rw
        # special case for the point at infinity
        if (exp == E(0)):
            ok = ok and res[2] == 0
        else:
            exp = (exp[0], exp[1])
            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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_mixed_add_projective_char2(E, a_2, a_6):
    """
    test mixed addition of points (affine and projective) on E2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    Fq = E.base_field()
    Fq1 = Fq(1)
    list_points = get_torsion_points(E, a_2, a_6)
    list_l = [Fq.fetch_int(item) for item in [1, 2, 3, 11]]
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        Rw = E(R[0], R[1]) # 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], S[1]) # long Weierstrass coordinates
                exp = Rw+Sw # long Weierstrass coordinates
                res = mixed_add_proj_char2(Rp, S, a_2, a_6)
                # special case for the point at infinity
                if (exp == E(0)):
                    ok = ok and res[2] == 0
                else:
                    exp = (exp[0], exp[1])
                    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], Rw[1], Fq1) # Montgomery coordinates
        S = (Sw[0], Sw[1], 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_char2(Rp, S, a_2, a_6) # Montgomery coordinates
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0], expw[1])
                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 + X3*Y3*Z3 == X3**3 + a_2 * X3**2 * Z3 + a_6 * Z3**3
                print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    print("test mixed_add_proj_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_dbl_projective_char2(E, a_2, a_6):
    """
    test doubling of a point in projective coordinates on E in characteristic 2
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    Fq = a_2.parent()
    list_points = get_torsion_points(E, a_2, a_6)
    list_l = [Fq.fetch_int(item) for item in [1, 2, 3, 11]]
    ok = True
    i = 0
    while ok and i < len(list_points):
        R = list_points[i]
        Rw = E(R[0], R[1]) # long Weierstrass coordinates
        expw = 2*Rw # 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
            res = dbl_proj_char2(Rp, a_2, a_6)
            # special case for the point at infinity
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0], expw[1])
                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 + X3*Y3*Z3 == X3**3 + a_2 * X3**2 * Z3 + a_6 * Z3**3
                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], Rw[1]) # 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_char2(Rp, a_2, a_6) # Montgomery coordinates
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0], expw[1])
                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 + X3*Y3*Z3 == X3**3 + a_2 * X3**2 * Z3 + a_6 * Z3**3
                print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    print("test dbl_proj_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_double_and_add_proj_char2(E, a_2, a_6):
    """
    test test_double_and_add_proj_char2(m, P, a_2, a_6)
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    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], Pw[1], Fq1)
        j = 0
        while ok and j < len(list_l):
            m = list_l[j]
            expw = m * Pw
            res = double_and_add_proj_char2(m, P, a_2, a_6)
            if (expw == E(0)):
                ok = ok and res[2] == 0
            else:
                exp = (expw[0], expw[1])
                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 + X3*Y3*Z3 == X3**3 + a_2 * X3**2 * Z3 + a_6 * Z3**3
                    print("output point on the curve: {}".format(point_on_curve))
            j = j+1
        i = i+1
    print("test double_and_add_proj_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_add_affine_x_only_char2(E, a_2, a_6):
    """
    test affine addition of points in x-only coordinates in characteristic 2
    It is assumed that the inputs are consistent, that is the points are distincts, non-opposite, and non-zero
    """
    assert E.a1() == 1 and E.a3() == 0 and E.a4() == 0 and E.a2() == a_2 and E.a6() == a_6
    # first, test addition of a torsion point with another distinct random point
    list_points = get_torsion_points(E, a_2, a_6)
    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], R[1], R[2])  # long Weierstrass coordinates
        S = (Sw[0], Sw[1], Sw[2]) # Montgomery coordinates
        Tw = Rw - Sw                  # long Weierstrass coordinates
        T = (Tw[0], Tw[1], Tw[2])
        x4 = T[0]
        x3 = add_affine_x_only_char2(R[0], S[0], x4)
        Qw = Rw + Sw
        Q = (Qw[0], Qw[1], 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], Rw[1], Fq1) # Montgomery coordinates
        S = (Sw[0], Sw[1], Fq1) # Montgomery coordinates
        T = (Tw[0], Tw[1], Fq1) # Montgomery coordinates
        x3 = add_affine_x_only_char2(R[0], S[0], T[0]) # Montgomery coordinates
        Qw = Rw+Sw # long Weierstrass coordinates
        Q = (Qw[0], Qw[1], 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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

def test_double_affine_x_only_char2(E, a_2, a_6):
    """
    test affine doubling of points in x-only coordinates in characteristic 2
    """
    list_points = get_torsion_points(E, a_2, a_6)
    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
        if R[0] == 0:
            i=i+1
            continue # it will skip the rest of the iteration and directly go to the next i
        Rw = E(R[0], R[1]) # long Weierstrass coordinates
        Qw = 2*Rw # long Weierstrass coordinates
        Q = (Qw[0], Qw[1], Qw[2])
        x3 = double_affine_x_only_char2(R[0], a_6)
        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], Rw[1], Fq1) # Montgomery coordinates
        x3 = double_affine_x_only_char2(R[0], a_6) # Montgomery coordinates
        Qw = 2*Rw
        Q = (Qw[0], Qw[1], Qw[2]) # Montgomery coordinates
        ok = x3 == Q[0]
        if not ok:
            print("Error double_affine_x_only_char2")
            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_char2 with a_6={}: {}".format(a_6.integer_representation(), ok))
    return ok

if __name__ == "__main__":
    # Define the base field (actually an extension of F2)
    m = 233
    R = PolynomialRing(GF(2), 'Z')
    Z = R.gen()
    F2m = GF(2**m, 'z', modulus = Z**233 + Z**74 + 1)

    # Define curve NIST B-233 with efficient halving
    a_2 = F2m(1)
    a_6 = F2m.fetch_int(0x0066647ede6c332c7f8c0923bb58213b333b20e9ce4281fe115f7d8f90ad)
    E2 = EllipticCurve(F2m,  [1, a_2, 0, 0, a_6])

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

    test_add_affine_char2(E2, a_2, a_6)
    test_dbl_affine_char2(E2, a_2, a_6)
    test_mixed_add_projective_char2(E2, a_2, a_6)
    test_dbl_projective_char2(E2, a_2, a_6)
    test_double_and_add_proj_char2(E2, a_2, a_6)

    test_add_affine_x_only_char2(E2, a_2, a_6)
    test_double_affine_x_only_char2(E2, a_2, a_6)
    test_add_proj_xz_only_char2(E2, a_2, a_6)
    test_mixed_add_proj_xz_only_char2(E2, a_2, a_6)
    test_double_proj_xz_only_char2(E2, a_2, a_6)
    test_montgomery_ladder_affine_x_only_char2(E2, a_2, a_6, r, cofactor)
    test_montgomery_ladder_proj_xz_only_char2(E2, a_2, a_6, r, cofactor)