def torsion_basis(E, l, e, p):
    n = E.order()
    cofactor = n // l^e
    while True:
        P = cofactor * E.random_point()
        if P * l^(e-1) != E(0):
            break
    while True:
        Q = cofactor * E.random_point()
        if Q * l^(e-1) != E(0):
            eP = P.weil_pairing(Q, l^e)
            if eP^(l^(e-1)) != eP.parent()(1):
                break
    return P, Q
def velu_isogeny(E, K):
    phi = E.isogeny(K)
    return phi, phi.codomain()

def find_curve_with_l_torsion(p, l):
    Fp = GF(p)
    for _ in range(2000):
        a, b = Fp.random_element(), Fp.random_element()
        if 4*a^3 + 27*b^2 == 0:
            continue
        E = EllipticCurve(Fp, [a, b])
        if E.order() % l == 0:
            return E
    raise ValueError(f"I cant find a curve over GF({p}) with {l}-torsion")

def random_l_isogeny(E, l):
    n = E.order()
    if n % l != 0:
        raise ValueError(f"l={l} does not divide #E = {n}")
    cofactor = n // l
    for _ in range(500):
        P = E.random_point()
        K = cofactor * P
        if K != E(0) and K.order() == l:
            return velu_isogeny(E, K)
    raise ValueError("no kernel point found after too many tries")


def csidh_params():
    small_primes = [3, 5, 7, 11, 13, 17, 19]
    product = 4 * prod(small_primes)
    for k in range(1, 10000):
        p = product * k - 1
        if p > 0 and is_prime(p) and p % 4 == 3:
            break
    else:
        raise ValueError("Could not find CSIDH prime")
    print(f"  [csidh_params] p = {p}  (p mod 4 = {p % 4})")
    Fp = GF(p)
    E0 = EllipticCurve(Fp, [1, 0])
    assert E0.is_supersingular(), "E0 must be supersingular"
    return p, small_primes, E0

def _csidh_step(E_cur, l, sign, p):
    Fp = GF(p)
    n = E_cur.order()
    if n % l != 0:
        E_try = E_cur.quadratic_twist()
        if E_try.order() % l == 0:
            E_cur = E_try
            n = E_cur.order()
        else:
            raise ValueError(f"neither E nor its twist has {l}-torsion (n={n}, are you fr now?)")
    cofactor = n // l
    for _ in range(1000):
        R = E_cur.random_point()
        K = cofactor * R
        if K != E_cur(0) and K.order() == l:
            phi = E_cur.isogeny(K)
            return phi.codomain()
    raise ValueError(f"no {l}-torsion kernel point found")

def csidh_action(E, exp_vec, primes, p):
    E_cur = E
    for l, e in zip(primes, exp_vec):
        if e == 0:
            continue
        sign = 1 if e > 0 else -1
        for _ in range(abs(e)):
            E_cur = _csidh_step(E_cur, l, sign, p)
    return E_cur

def make_sidh_prime(eA, eB):
    for f in range(1, 100000):
        p = 2^eA * 3^eB * f - 1
        if p > 0 and is_prime(p):
            return p, f
    raise ValueError("I cant find a SIDH prime!!!!!!!!!!")

def sidh_setup(p, eA, eB):
    R.<x> = GF(p)[]
    Fp2.<i> = GF(p^2, modulus=x^2+1)
    E0 = EllipticCurve(Fp2, [1, 0])
    PA, QA = torsion_basis(E0, 2, eA, p)
    PB, QB = torsion_basis(E0, 3, eB, p)
    return E0, PA, QA, PB, QB

def sidh_keygen_A(E0, PA, QA, PB, QB, eA):
    while True:
        mA, nA = randint(0, 2^eA-1), randint(0, 2^eA-1)
        if (mA, nA) != (0, 0):
            break
    KA = mA*PA + nA*QA
    phi_A = E0.isogeny(KA, algorithm='factored')
    EA = phi_A.codomain()
    return (mA, nA), EA, phi_A(PB), phi_A(QB)

def sidh_keygen_B(E0, PA, QA, PB, QB, eB):
    while True:
        mB, nB = randint(0, 3^eB-1), randint(0, 3^eB-1)
        if (mB, nB) != (0, 0):
            break
    KB = mB*PB + nB*QB
    phi_B = E0.isogeny(KB, algorithm='factored')
    EB = phi_B.codomain()
    return (mB, nB), EB, phi_B(PA), phi_B(QA)

def sidh_shared_A(mA, nA, EB, phi_PA, phi_QA):
    KA = mA*phi_PA + nA*phi_QA
    return EB.isogeny(KA, algorithm='factored').codomain().j_invariant()

def sidh_shared_B(mB, nB, EA, phi_PB, phi_QB):
    KB = mB*phi_PB + nB*phi_QB
    return EA.isogeny(KB, algorithm='factored').codomain().j_invariant()
def isogeny_graph_neighbors(E, l):
    poly = E.division_polynomial(l)
    neighbors = []
    seen_j = set()
    for fac, _ in poly.factor():
        if fac.degree() != 1:
            continue
        x0 = -fac[0] / fac[1]
        try:
            pts = E.lift_x(x0, all=True)
            for P in pts:
                if P != E(0) and (l * P) == E(0):
                    phi = E.isogeny(P)
                    F = phi.codomain()
                    j = F.j_invariant()
                    if j not in seen_j:
                        seen_j.add(j)
                        neighbors.append(F)
                    break
        except ValueError:
            pass
    return neighbors

def bfs_isogeny_path(E_start, j_end, l, max_depth=8):
    from collections import deque
    queue = deque([(E_start, [E_start.j_invariant()])])
    visited = {E_start.j_invariant()}
    while queue:
        E, path = queue.popleft()
        if E.j_invariant() == j_end:
            return path
        if len(path) > max_depth:
            continue
        for F in isogeny_graph_neighbors(E, l):
            j = F.j_invariant()
            if j not in visited:
                visited.add(j)
                queue.append((F, path + [j]))
    return None