Velu's formula for computing isogeny in short Weirstrass form

velu.sage

from sage.all_cmdline import *

Fq=GF(11)
E=EllipticCurve(Fq, [3,7])
G0=(E.gens()[0])
R = PolynomialRing(Fq, 'X,Y'); R.inject_variables()

def normalize_rational_map(rational_map, curve : EllipticCurve):
  def reduce_y2(poly, fx):
    xterms = 0
    yterms = 0
    for (i,xcoeff) in enumerate(poly.polynomial(Y).coefficients(sparse=False)):
      if i % 2 == 0:
        xterms = xterms + xcoeff*(fx ** (i//2))
      else:
        yterms = yterms + xcoeff*(fx ** ((i-1)//2))
    return (xterms , yterms)

  fx = X**3 + curve.a4()*X + curve.a6();
  (n_x,n_y) = reduce_y2(rational_map.numerator(), fx)
  (d_x,d_y) = reduce_y2(rational_map.denominator(), fx)

  if d_y != 0:
    den = d_x*d_x - (d_y*d_y*fx)
    # num = (n_x + Y*n_y)*(d_x - Y*d_y)
    num = (n_x*d_x - n_y*d_y*fx) + Y*(n_y*d_x - n_x*d_y)
  else:
    den = d_x
    num = n_x + Y*n_y

  return (num/den)

def translation_map(P):
  m = (Y-P[1])/(X-P[0])
  x_t = m*m - P[0] - X
  y_t = m*(P[0] - x_t) - P[1]
  return (normalize_rational_map(x_t, P.curve()), 
          normalize_rational_map(y_t, P.curve()))

def phi(points):
  xs = X
  ys = Y
  curve = None
  for p in points:
    if p[0] == 0 and p[1] == 1 and p[2] == 0:
      continue
    curve = p.curve()
    xp,yp = translation_map(p)
    xs = xs + (xp - p[0])
    ys = ys + (yp - p[1])
  return (normalize_rational_map(xs,curve),normalize_rational_map(ys,curve))


(phi_x,phi_y)=phi([i*G0 for i in range(1,G0.order())])
print(phi_x.denominator().factor())
print([phi_x.denominator()((i*G0)[0],(i*G0)[1])  for i in range(1,G0.order())])