Q := RationalField();
JSeq := [0, 54000, -12288000, 1728, 287496, -3375, 16581375, 8000, -32768, 
-884736, -884736000, -147197952000, -262537412640768000];

F<b,c> := FunctionField(Q,2);

E := EllipticCurve([F | 1-c, -b, -b, 0, 0]);
EP := E![0,0];

j := jInvariant(E);
Nj := Numerator(F!j);
Dj := Denominator(F!j);

P<x,y> := PolynomialRing(Q,2);

PJSeq := [P!Nj - j*P!Dj : j in JSeq];

/* These are the 13 polynomials we get when we set j(b,c) = j_i, for each 
of the j_i in the sequence JSeq. They're written as polynomials in x and 
y. */




TC := function(i,N,upplim) 

/* This function gives degree sequences mod primes from 3 to upplim of the 
factorization of the one-variable polynomial whose roots are the c-values 
of elliptic curves with an N-torsion point. For now I'm taking N to be 
odd. */

x1 := (((N-1) div 2)*EP)[1];
x2 := (((N+1) div 2)*EP)[1];

/* These are the x-coordinates of plus or minus (N-1)/2 times the origin, 
and plus or minus (N+1)/2 times the origin. If they coincide, then a 
little thought shows that N times the origin is the point at infinity. */

Np := P!(Numerator(x1)*Denominator(x2) - Numerator(x2)*Denominator(x1));

Ri := Resultant(PJSeq[i], Np, y);


Gi := Lcm([Denominator(c) : c in Coefficients(Ri)]);
Ri := Ri*Gi;
Li := Gcd([Numerator(c) : c in Coefficients(Ri)]);
Ri := Ri*Li;



/* I suppose that the way to speed the process up is to precompute the 
resultant, before picking which j-invariant we use, and then to define Ri 
to be the general resultant with j_i plugged in. */

p := 3;
mindeg := 1;

  while p le upplim do

  FF := Factorization(PolynomialRing(GF(p))!UnivariatePolynomial(Ri));
  degseq := [Degree(FF[i][1]) : i in [1 .. #FF]];

  if degseq[1] gt mindeg then
  mindeg := degseq[1];
  end if;
  
  p, degseq;

  p := NextPrime(p);
  end while;

return "Looks like the smallest degree is", mindeg;

end function;


// So for example you might try running TC(1,31,200).