_<t> := FunctionField(Rationals());
QQt := Parent(Numerator(t));
_<x> := PolynomialRing(Parent(t));

// Tate form <ns, b, c, Ps>
dat := [<[3], 0, 0, [[]]>,
        <[4], t, 0, [[0,0]]>,
        <[2,2], 0, 0, [[]]>,
        <[5], t, t, [[0,0]]>,
        <[6], t^2+t, t, [[0,0]]>,
        <[7], t^3-t^2, t^2-t, [[0,0]]>,
        <[8], 2*t^2 - 3*t + 1, (2*t^2 - 3*t + 1)/t, [[0,0]]>,
        <[2,4], (t^2-1)/16, 0, [[(t-1)/8, (t-1)^2/32], [0,0]]>,
        <[9], t^5 - 2*t^4 + 2*t^3 - t^2, t^3 - t^2, [[0,0]]>,
        <[10], (2*t^5 - 3*t^4 + t^3)/(t^4 - 6*t^3 + 11*t^2 - 6*t + 1), 
             (-2*t^3 + 3*t^2 - t)/(t^2 - 3*t + 1), [[0,0]]>,
        <[12], (12*t^6 - 30*t^5 + 34*t^4 - 21*t^3 + 7*t^2 - t)/
             (t^4 - 4*t^3 + 6*t^2 - 4*t + 1),
             (-6*t^4 + 9*t^3 - 5*t^2 + t)/(t^3 - 3*t^2 + 3*t - 1),
             [[0,0]]>,
        <[2,6], (1-t^2)/(6*t-10)+((1-t^2)/(6*t-10))^2, (1-t^2)/(6*t-10), 
             [[(1/18*t^3 - 7/18*t^2 + 5/6*t - 1/2)/(t^2 - 10/3*t + 25/9),
               (1/108*t^5 - 11/108*t^4 + 23/54*t^3 - 5/6*t^2 + 3/4*t - 1/4)/(t^3 - 5*t^2 + 25/3*t - 125/27)], [0,0]]>, 
        <[2,8], (t^4 - 28*t^3 + 304*t^2 - 1472*t + 2560)/(t^4 - 16*t^3 + 32*t^2 + 256*t + 256),
        (1/4*t^4 - 7*t^3 + 76*t^2 - 368*t + 640)/(t^3 - 20*t^2 + 80*t + 192),
        [[(-1/64*t^4 + 1/2*t^3 - 6*t^2 + 32*t - 60)/(t^2 - 8*t - 16),
          (-1/512*t^7 + 13/128*t^6 - 73/32*t^5 + 229/8*t^4 - 431/2*t^3 + 966*t^2 - 2360*t + 2400)/
          (t^4 - 16*t^3 + 32*t^2 + 256*t + 256)], [0,0]]>];

for i := 1 to #dat do
  ns, b, c, Ptates := Explode(dat[i]);

  if ns eq [3] then
    Etate := EllipticCurve([1,0,t,0,0]);
    Ptates := [Etate![0,0]];
  elif ns eq [2,2] then
    Etate := EllipticCurve(x*(x-1)*(x-t));
    Ptates := [Etate![0,0], Etate![1,0]];
  else
    Etate := EllipticCurve([1-c,-b,-b,0,0]);
    Ptates := [Etate!P : P in Ptates];
  end if;
  assert &and[ns[j]*Ptates[j] eq Etate!0 and
                 &and[d*Ptates[j] ne Etate!0 : d in Divisors(ns[j]) | d ne ns[j]] : j in [1..#Ptates]];

  Ew0, mEw0 := WeierstrassModel(Etate);
  dfull := Denominator(Coefficients(Ew0)[4]);
  bl, d := IsPower(dfull, 4);
  assert bl and Denominator(Coefficients(Ew0)[5]) eq d^6;
  E := QuadraticTwist(Ew0, d^2);
  Ps := [E![mEw0(Ptates[j])[1]*d^2, mEw0(Ptates[j])[2]*d^3] : j in [1..#Ptates]];
  assert &and[ns[j]*Ps[j] eq E!0 : j in [1..#Ps]];

  G := AbelianGroup(ns);
  for Hrec in Subgroups(G) do
    H := Hrec`subgroup;
    if #H eq 1 then 
      Ep := E;
    else
      HPs := [ &+[Eltseq(G!h)[j]*Ps[j] : j in [1..#Eltseq(G!h)] ] : h in H];
      Ep, mp := IsogenyFromKernel(E, &*[x-HP[1] : HP in HPs | HP ne E!0]);
    end if;
    printf "f = %o, g = %o", Coefficients(Ep)[4], Coefficients(Ep)[5];
    rs := [Degree(QQt!c) : c in Coefficients(Ep)[4..5]];
    noverm := Max(rs[1]/4,rs[2]/6);
    assert Denominator(noverm) eq 1 or Numerator(noverm) eq 1;
    powerofX := (Denominator(noverm)+1)/(12*Numerator(noverm));
    printf "HS: For %o and subgroup %o, leadterm %o, r = %o, s = %o => X^(%o)",
          ns, Invariants(H), 
          LeadingTerm(QQt!Coefficients(Ep)[4]) eq LeadingTerm(QQt!Coefficients(E)[4]) and
          LeadingTerm(QQt!Coefficients(Ep)[5]) eq LeadingTerm(QQt!Coefficients(E)[5]),
          rs[1], rs[2], powerofX;

    // same, but with quadratic twists
    noverm := Max(rs[1]/2,rs[2]/3);
    if Denominator(noverm) eq 1 or Numerator(noverm) eq 1 then
      n := Numerator(noverm);
      m := Denominator(noverm);
      if m+1 gt n then
        powerofXquadtwist := (Denominator(noverm)+1)/(6*Numerator(noverm));      
      elif m+1 eq n then
        powerofXquadtwist := "1/6+eps";      
      else
        powerofXquadtwist := 1/6;
      end if;
      assert Gcd([QQt | c : c in Coefficients(Ep)]) eq 1;
      printf ", or X^(%o) if quadtwist\n", powerofXquadtwist;
    else
      printf ", no quad twist?\n";
    end if;
  end for;
end for;

