A<x,y> := PolynomialRing(GF(2),2);
P,homogenize := Homogenization(A);
print("f := The polynomial from math.mit.edu/~drew/FFFc53.txt reduced mod 2");
f := y^12*(y+1)^4*x^39+y^12*(y+1)^5*x^38+y^6*(y^18+y^15+y^11+y^10+y^5+y^4+y^3+y^2+1)*(y^4+y^3+y^2+y+1)*(y+1)^3*x^25+y^5*(y^10+y^9+y^8+y^6+y^2+y+1)*(y^15+y^14+y^12+y^10+y^9+y^5+y^2+y+1)*(y+1)^3*x^23+(y^3+y^2+1)*(y^23+y^22+y^20+y^18+y^16+y^15+y^9+y^8+y^5+y^4+1)*y^5*(y+1)^3*x^22+(y^3+y^2+1)^2*y^11*(y^4+y^3+1)*(y^2+y+1)*x^33+(y^3+y^2+1)^2*(y^8+y^7+y^5+y+1)*y^8*(y+1)^3*x^31+y^7*(y^5+y^4+y^3+y^2+1)*(y^11+y^9+y^7+y^4+1)*(y+1)^5*x^28+y^8*(y^5+y^3+y^2+y+1)*(y^7+y^6+y^5+y^4+y^2+y+1)*(y+1)^9*x^27+y^11*(y^6+y^5+y^4+y+1)*x^36+y^5*(y^31+y^29+y^28+y^25+y^22+y^21+y^20+y^19+y^18+y^17+y^16+y^13+y^12+y^6+y^2+y+1)*x^20+(y^11+y^9+y^7+y^6+y^4+y^2+1)*(y+1)^14*x^5+y^5*(y+1)^16*x+(y^5+y^3+1)*y^8*(y^2+y+1)*(y^7+y^6+y^4+y^2+1)*(y+1)^8*x^26+y*(y^5+y^4+y^3+y^2+1)^2*(y^13+y^7+y^5+y^3+y^2+y+1)*(y+1)*(y^9+y^8+y^7+y^3+y^2+y+1)*x^16+y*(y^22+y^21+y^20+y^19+y^18+y^17+y^11+y^10+y^9+y^8+y^5+y^2+1)*(y^2+y+1)*(y^6+y^3+1)*(y+1)^4*x^15+(y^3+y+1)*y*(y^17+y^16+y^14+y^12+y^7+y^5+y^4+y+1)*(y^7+y^6+y^5+y^4+1)*(y+1)^5*x^14+y*(y^15+y^14+y^11+y^6+1)*(y+1)^6*(y^5+y^4+y^3+y+1)*(y^6+y^5+y^4+y^2+1)*x^13+y*(y^5+y^2+1)*(y^7+y^5+y^4+y^3+y^2+y+1)*(y+1)^9*(y^14+y^9+y^8+y^7+y^5+y^3+y^2+y+1)*x^10+y^5*(y^2+y+1)*(y+1)^22*x^3+(y^15+y^14+y^13+y^12+y^10+y^8+y^7+y^5+y^4+y^2+1)*y^8*(y+1)^3*x^30+(y^8+y^6+y^5+y^4+y^2+y+1)*y^7*(y+1)^12*x^29+y^5*(y^24+y^23+y^22+y^21+y^19+y^18+y^16+y^15+y^12+y^11+y^10+y^9+y^8+y^6+y^5+y^4+1)*(y+1)^3*x^24+(y^18+y^16+y^15+y^14+y^12+y^11+y^8+y^5+y^2+y+1)*y^4*(y^14+y^13+y^12+y^11+y^10+y^9+y^8+y^6+y^5+y^4+1)*x^18+(y^25+y^24+y^23+y^22+y^19+y^18+y^17+y^16+y^15+y^14+y^12+y^11+y^4+y+1)*y^2*(y^8+y^7+y^6+y^4+y^2+y+1)*x^17+y^11*(y^2+y+1)*(y+1)^4*x^37+(y^3+y^2+1)*y^12*(y^5+y^4+y^3+y^2+1)*x^35+y^12*(y^5+y^4+y^3+y^2+1)*(y^5+y^4+y^2+y+1)*x^34+y^11*(y^12+y^11+y^8+y^2+1)*(y+1)*x^32+y^4*(y+1)^17+y^5*(y^5+y^4+y^3+y^2+1)*(y^21+y^20+y^17+y^15+y^14+y^12+y^9+y^8+y^7+y^6+y^4+y^2+1)*(y+1)^4*x^21+y^4*(y^15+y^12+y^11+y^8+y^6+y^5+y^4+y^3+y^2+y+1)*(y^6+y^3+1)*(y^12+y^11+y^10+y^9+y^7+y^5+y^4+y^3+1)*x^19+(y^12+y^11+y^10+y^7+y^6+y^4+y^2+y+1)*y*(y^16+y^14+y^9+y^6+y^3+y+1)*(y+1)^7*x^12+(y^12+y^11+y^8+y^7+y^6+y^5+y^4+y^3+1)*y*(y^16+y^13+y^12+y^11+y^10+y^8+y^5+y+1)*(y+1)^8*x^11+(y^5+y^3+1)*(y^13+y^12+y^11+y^9+y^7+y^5+y^4+y^3+1)*y*(y+1)^10*x^9+(y^3+y^2+1)*(y^14+y^12+y^10+y^8+y^6+y^5+y^4+y^2+1)*y*(y+1)^11*x^8+(y^14+y^13+y^11+y^8+y^6+y^5+y^4+y^2+1)*y*(y^2+y+1)*(y+1)^12*x^7+y^2*(y^6+y^4+y^3+y+1)*(y^4+y+1)*(y+1)^13*x^6+(y^3+y^2+1)*y^4*(y^6+y+1)*(y+1)^15*x^4+(y^3+y^2+1)*(y^3+y+1)*y^6*(y+1)^16*x^2 ;

f := homogenize(f);
XF2 := Scheme(ProjectiveSpace(P),f);
FF := FunctionField(XF2);
AFF := AlgorithmicFunctionField(FF);
plc1 := Places(AFF,1);
cuspsum := &+ plc1;

print("This file will show that X1(53)(F2) has no functions of degree <= 16");
print("so the F2-gonality and hence the Q-gonality of X1(53) is at least 17");

print("The number of places of X1(53)(F2) of degrees 1..7 is:");
[#Places(AFF,i) : i in [1..7]];

// Result: [ 26, 0, 0, 0, 0, 0, 0 ]

// Say f is a function of degree <= 16, since f^(-1)({0,1,infinity}) contains
// all 26 rational points, by the pigeon hole principle we may assume that
// f^(-1)({infinity} contains at least 9 distinct rational points, and so
// the remaining polar-part of f has degree at most 7.

// Either there is one multiplicity > 3, which is covered by this:

Dimension(RiemannRochSpace( 3*cuspsum + 5*plc1[1] ));

// or there are two.  Then one of them is 3+2 and the other is 3+1.
// We fix the +2 at plc1[1].  The other has to loop through plc1[2..26] and
// we do this covering 5 at the same time:

{Dimension(RiemannRochSpace( 3*cuspsum + 2*plc1[1] + &+plc1[ 5*i - 3 .. 5*i + 1 ]  )) : i in [1..5]};
