#### computes isogeny following the BT directions

def compute_isogeny_BT_direction(EC, ell, i, j, points):
    """
    computes the ell-isogeny from the curve E in the BT direction i, 
    using the points P,Q (images of the basis of the starting curve) and T (the previous kernel approximation)
    """

    P, Q, T = points

    if i<  ell:
        S = (i*ell^(j))*P + T
        n = S.order()
    
    else:
        S = Q
        n = S.order()
    
    phi = EC.isogeny(int(n/ell)*S)
    
    new_points = [phi(P), phi(Q), phi(S)]
    return phi, new_points


def torsion_basis_large_N(E, N):
    """
    outputs any basis of E[N]
    by finding independent random points of order N
    """
    found = False
    
    while not found:
        xP= R.random_element()
        try:
            P = E.lift_x(xP)
            P = int((p+1)/N)*P
            if P.order() == N:
                found = True
        except:
            continue

    found = False
    while not found:
        xQ= R.random_element()
        try:
            Q = E.lift_x(xQ)
            Q = int((p+1)/N)*Q
            if (Q.order() == N and P.weil_pairing(Q, N)^(ell^(k-1)) != 1):
                found = True
        except:
            continue
            
    return P,Q


def path_to_isogeny(E, ell, basis, path):
    """
    given a path in the BT tree as a list of strings (i1, ..., ik)
    and a basis P, Q of E[ell^k]
    constructs the corresponding path of elliptic curves E0 = E, E1, ..., Ek
    
    Assume i1 != infty:
    By construction, Ek = E/< sum ij*ell^(j-1) * P+Q>
    """
    
    path = [a for a in steps]
    path.reverse()
    print path
    
    curve = E
    P,Q = basis
    S = Q
    points = [P,Q,S]
    
    for j in range(len(path)):
        i = path.pop()
        phi, new_points = compute_isogeny_BT_direction(curve, ell, i, j, points)
        curve = phi.codomain()
        points = new_points
        print "direction", i, "curve", curve.j_invariant()
        

########################      
#### choose your ell
ell = 2


################## 
### set up a prime 3 mod 4 that has sufficiently large ell-torsion to be able to walk far

p = next_prime(500)
while (p % 4 != 3 or valuation(p+1, ell) < 8):
    p = next_prime(p)
    
print 'prime', p 

k = valuation(p+1, ell)

print 'valuation', k

# fixing the finite field once and for all
R.<a> = GF(p^2 )
E = EllipticCurve_from_j(287496).change_ring(R)
print 'starting elliptic curve', E

### find an independent basis
### comment out to repeat the walks

basis= torsion_basis_large_N(E, ell^k)
P,Q = basis

#### sequence of steps
steps = [1,0,0,0]

print steps

path_to_isogeny(E, ell,  basis, steps)


### outputs the j-invariant of the final curve
def check(E, ell,  basis, steps):
    P, Q = basis
    T = Q
    for j in range(len(steps)):
        T =T + steps[j]*ell^j*P
        
    if len(steps) < k:
        T = (ell^(k-len(steps)))*T
    Etarget = E.isogeny(  T).codomain()
    return Etarget.j_invariant()


check(E, ell,  basis, steps)