Find basis for the set of torsion points E[m]

Question

In paper "On the Cost of Computing Isogenies Between Supersingular Elliptic Curves" (source) reads

Let ${P, Q}$ be a basis for $E[2^{e/2}]$. Let $R_0 = [2^{e/2}−1]P , R_1 = [2^{e/2}−1]Q, R_2 = R_0 + R_1$ be the order-2 points on E.

This involves randomizing points, multiplying them by a factor $(p+1)/2^{e/2}$ as well as checking if their Weil-pairing is an $e/2$-th root of unity (see link1 and Chapman's answer).

Actually it is not a trivial task to find such basis on an arbitrary curve $E$ such that $R_0, R_1$ and $R_2$ cover all three points of order $2$ on $E$. Practically, computing $R_0, R_1$ and $R_2$ from randomized $P$ and $Q$ often results in at most two distinct points of order 2.

Does someone know a more efficient way to generate the basis points $P$ and $Q$ in this case? Thanks a lot!

Updated: Below is the algorithm in the paper Efficient compression of SIDH public keys (Sec. 3.2) implemented in SageMath

# this parameter set not work!
eA=8
eB=5
cst=1

# but this works!
#eA=4
#eB=3
#cst=1

# characteristic p in form c*2^eA*3^eB-1
p=cst*2**eA*3**eB-1 
_.<I> = GF(p)[]
K.<i> = GF(p^2, modulus=I^2+1)
E = EllipticCurve(K, [0, 6, 0, 1, 0])  # define base curve E

def find_torsion_basis(basic_curve, exponent,A,j):
    u=4+i
    global cst
    while True:
        x_R=j*u
        y_R_sqr=x_R**3 + A*x_R**2 + x_R

        if y_R_sqr.is_square():
            y_R=sqrt(y_R_sqr)
            basis_point=basic_curve([x_R,y_R])
            basis_point=cst*3**eB*basis_point
            break
        else:
            j=j+1
    return (j,basis_point)

A=6
(j,P_e)=find_torsion_basis(E,eA,A,1)
(j,Q_e)=find_torsion_basis(E,eA,A,j+1)
assert(P_e.weil_pairing(Q_e, Integer(2**(eA))).multiplicative_order()==2**(eA))

As an example, it works for $e_A = 4, e_B=3, c=1$ but NOT for $e_A = 8, e_B=5, c=1$, given $p$ of the form $p = c \cdot 2^{e_A} \cdot 3^{e_B} - 1$. For $p = 2^8 \cdot 3^5 -1$, the outputs $P$ and $Q$ don't even satisfy the Weil-pairing condition for being torsion basis of $E[2^{e_A}]$.

Please correct if I'm doing something wrong!

Updated2:

An important point as @BenSmith remarked is $u$ must be initialized as a non-square element $\mathbb{F}_{p^2}$. Thus, I adapted my code such that $u$ is the $x$-coordinate of a random point in $E(\mathbb{F}_{p^2})$ and tested for the condition that $u$ must be a non-square as mentioned in Efficient compression of SIDH public keys (Sec. 3.1).

def find_torsion_basis(basic_curve,u,A,j):
    global cst
    
    while True:
        x_R=j*u
        y_R_sqr=x_R**3 + A*x_R**2 + x_R

        if y_R_sqr.is_square():
            y_R=sqrt(y_R_sqr)
            basis_point=basic_curve([x_R,y_R])
            basis_point=cst*3**eB*basis_point
            break
        else:
            j=j+1
    return (j,basis_point)

def extract_real_and_img_from_complex_num(c): # return (a,b) given (a+b*i)
    c=str(c)
    components=[s for s in c.split("+")]
    for x in components:
        if "i" in x:
            img = int(x[:x.index("i")-1])
        else:
            real = int(x)
    return (real,img)

A=6
#find u, which is non-square in F_p^2
while True:
    u=E.random_point()[0]
    (real,img)=extract_real_and_img_from_complex_num(u)
    tmp=R(real)**2+ R(img)**2
    if not tmp.is_square(): # if tmp is non-square in F_p
        break
assert(not u.is_square())       
print("u =",u)

(j,P_e)=find_torsion_basis(E,u,A,1)
(j,Q_e)=find_torsion_basis(E,u,A,j+1)
assert(P_e.weil_pairing(Q_e, Integer(2**(eA))).multiplicative_order()==2**(eA))

Eventually it works! But sometimes the code above, given $u$ non-square (sure because the first assert never fires), still produces a wrong pair $P_e$ and $Q_e$, which triggers the assert on the last line.

Answer 1

Generating 2-power torsion bases (deterministically) is an important step in SIDH key compression. The idea is to construct a 2-torsion basis - that is, start from $R_0$ and $R_1$ - and then find points $P$ and $Q$ with nontrivial 2-Tate pairing with $R_0$ and $R_1$, respectively. The 2-Tate pairing has a nice explicit form which makes this easy to do. The general details, and an efficient method for doing this (especially for 2-power torsion bases), are given in Section 3 of Costello, Jao, Longa, Naehrig, Renes, and Urbanik's Efficient compression of SIDH public keys.

Edit: Costello et al. gives a simple method of generating points of full 2-power order, but it doesn't guarantee that two points generated this way will form a basis - you still need to check the Weil pairing. But just as Costello et al. forces nontrivial Tate pairings with $(0,0)$, you can go one step further and use the Tate pairing with the other two nontrivial two-torsion points to simultaneously construct two points generating $E(\mathbb{F}_{p^2})/[2]E(\mathbb{F}_{p^2})$. Multiplying these two ponts by the cofactor will yield a basis of the 2-power torsion. This is done in Section 3 of Zanon, Simplicio, Pereira, Doliskani, and Barreto's Faster key compression for isogeny-based cryptosystems.