nTorsonIsogeny

def normalize(polyList, fx):
    const_term  = 0
    y_term = 0

    for i, v in enumerate(polyList):
        xterm = polyList[i]
        if i % 2 == 0:
            expo = int(i/2)
            const_term = const_term + xterm * (fx ^ expo) 
        else:
            expo = int(int(i-1)/2)
            y_term = y_term + xterm * (fx ^ expo)

    return (y_term,const_term)


def nTorsionIsogney(curve, c):
    A,B          = curve.a4(), curve.a6()
    fp           = curve.base_ring()
    divP_Plus_2  = curve.division_polynomial(c+2,two_torsion_multiplicity=1)
    divP_Plus_1  = curve.division_polynomial(c+1,two_torsion_multiplicity=1)
    divP_minus_1 = curve.division_polynomial(c-1,two_torsion_multiplicity=1)
    divP_minus_2 = curve.division_polynomial(c-2,two_torsion_multiplicity=1)
    divP         = curve.division_polynomial(c,two_torsion_multiplicity=1)
    
    if c % 2 == 0:
        (X,Y)        = divP.variables()
    else:
        (X,Y)        = divP_Plus_1.variables()
        
    numPhi       = (X*divP*divP - divP_Plus_1*divP_minus_1).polynomial(Y).list()
    denPhi       = (divP*divP).polynomial(Y).list()
    numOmega     = (divP_Plus_2*divP_minus_1*divP_minus_1 - 
                    divP_minus_2*divP_Plus_1*divP_Plus_1).polynomial(Y).list()
    denOmega     =  (4*Y*divP*divP*divP).polynomial(Y).list()
    
    fx           = X^3+A*X+B;
    (numX,numCX)  = normalize(numPhi, fx)
    (denX,denCX)  = normalize(denPhi, fx)
    (numY,numCY)  = normalize(numOmega, fx)
    (denY,denCY)  = normalize(denOmega, fx)
    return (((numCX + Y*numX),(denCX + Y*denX)), 
            ((numCY + Y*numY),(fx*(denCY + Y*denY)/Y)))

Fp = GF(43)
FpBar=Fp.algebraic_closure();
E  = EllipticCurve(Fp, [1,3])
g = E.gens()[0];
kTor = 8;
((numX,denX), (numY, denY))=nTorsionIsogney(E, kTor)
numX=numX.change_ring(FpBar)
print (numX(g[0],g[1])/denX(g[0],g[1]), (g[1]*numY(g[0],g[1])/denY(g[0],g[1])))