import sys

print('version 2020.03.17')

for q in [7,9,11,13,17,19,23,25,27,29,31]:
  k.<fieldgen> = GF(q)
  kX.<X> = k[]

  for A in k:
    a1 = 0
    a2 = A
    a3 = 0
    a4 = 1
    a6 = 0
    try:
      E = EllipticCurve([a1,a2,a3,a4,a6])
      if E.cardinality().is_power_of(2): continue
    except:
      continue

    for P in E:
      l = P.order()
      if not l.is_prime(): continue
      if l == 2: continue

      print('try','isogeny',q,A,P,l)
      sys.stdout.flush()

      x = {s:(s*P)[0] for s in range(1,l)}
      S = range(1,l-1,2)
      hS = kX(prod(X-x[s] for s in S))

      # renes formulas for codomain curve:

      c0 = prod(x[s] for s in range(1,(l-1)/2+1))
      assert c0 == prod(x[s] for s in S)
      assert c0^2 == prod(x[s] for s in range(1,l))
      assert c0 == (-1)^((l-1)//2)*hS(0)
      sigma = sum(x[s]-1/x[s] for s in range(1,l))
      newA = c0^2*(A-3*sigma)

      # moody-shumow formulas for newA via edwards:

      r = hS(1)/hS(-1)
      assert r == prod((x[s]-1)/(x[s]+1) for s in S)
      d = ((A-2)/(A+2))^l * r^8
      assert newA == 2*(1+d)/(1-d)

      # renes formulas for newX, newY:

      newX = X*prod((x[s]*X-1)/(X-x[s]) for s in range(1,l))
      assert newX == X*prod((x[s]*X-1)/(X-x[s]) for s in S)^2
      assert newX == X^l*hS(1/X)^2/hS(X)^2

      newXnumer = X*prod(x[s]*X-1 for s in range(1,l))
      newXdenom = prod(X-x[s] for s in range(1,l))
      assert newX == newXnumer/newXdenom
      newXdiff = (newXdenom*newXnumer.diff()-newXnumer*newXdenom.diff())/newXdenom^2

      YY = X^3+A*X^2+X
      newYY = c0^2*YY*newXdiff^2
      assert newYY == newX^3+newA*newX^2+newX

      # newA via generic point of order 2:

      R = kX.quotient(X^2+A*X+1)
      alpha = R(X)
      newalpha = alpha^l*hS(1/alpha)^2/hS(alpha)^2
      assert newA == -(newalpha+1/newalpha)

      # newX, newY via jets:

      kjets.<epsilon> = k[[]]
      kjetsY.<Y> = kjets[]

      for alpha in k:
        if alpha == 0: continue
        if hS(alpha) == 0: continue

        alphajet = alpha+epsilon+O(epsilon^2)

        newalphajet = alphajet^l*hS(1/alphajet)^2/hS(alphajet)^2
        assert newalphajet == alphajet*hS.reverse()(alphajet)^2/hS(alphajet)^2
        newalpha = newalphajet[0]
        newalphadiff = newalphajet[1]

        assert newalpha == newX(alpha)
        assert newalphadiff == newXdiff(alpha)

        beta2 = alpha^3+A*alpha^2+alpha
        newbeta2 = c0^2*beta2*newalphadiff^2
        assert newbeta2 == newalpha^3+newA*newalpha^2+newalpha

      # fast evaluation:

      def F0(X1,X2):
        return (X1-X2)^2

      def F1(X1,X2):
        return -2*((X1*X2+1)*(X1+X2)+2*A*X1*X2)

      def F2(X1,X2):
        return (X1*X2-1)^2

      b = floor(sqrt(l-1)/2)
      if b == 0:
        b2 = 0
      else:
        b2 = floor((l-1)/(4*b))
      I = set(2*b*(2*i+1) for i in range(b2))
      J = set(range(1,2*b,2))
      K = range(4*b*b2+1,l-1,2)

      for alpha in k:
        hI = kX(prod(X-x[i] for i in I))
        EJ = kX(prod(F0(X,x[j])*alpha^2+F1(X,x[j])*alpha+F2(X,x[j]) for j in J))
        denom = hI.resultant(EJ)*kX(prod(alpha-x[s] for s in K))

        if denom == 0:
          assert alpha != 0
          assert hS(alpha) == 0
          continue

        if alpha == 0:
          assert newX(alpha) == 0
          continue

        EJ2 = kX(prod(F0(X,x[j])/alpha^2+F1(X,x[j])/alpha+F2(X,x[j]) for j in J))
        numer = hI.resultant(EJ2)*kX(prod(1/alpha-x[s] for s in K))

        assert alpha^l*(numer/denom)^2 == newX(alpha)