Q := RationalField(); JSeq := [0, 54000, -12288000, 1728, 287496, -3375, 16581375, 8000, -32768, -884736, -884736000, -147197952000, -262537412640768000]; F := 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 := 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 := 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).