A<x,y> := PolynomialRing(GF(2),2);
P,homogenize := Homogenization(A);
print("f := An improved equation for X1(57) reduced mod 2");
f := (y^2+y+1)*(y^3+y+1)*(y+1)^2*(y^6+y^4+y^2+y+1)*y^5*(y^5+y^3+1)*x^19+(y^2+y+1)^15*x+(y+1)^8*x^26+(y+1)^10*x^30+(y^2+y+1)^15+(y^2+y+1)*(y+1)^6*(y^8+y^4+y^3+y^2+1)*x^27+(y+1)^8*(y^8+y^5+y^3+y+1)*(y^5+y^4+y^2+y+1)*x^24+(y^16+y^14+y^11+y^10+y^9+y^8+y^7+y^6+y^5+y+1)*(y^2+y+1)*(y+1)^5*x^21+(y^16+y^11+y^10+y^9+y^7+y^4+1)*(y+1)^3*(y^3+y^2+1)*x^20+(y^2+y+1)^2*(y^12+y^11+y^10+y^8+1)^2*y^2*x^14+(y^13+y^11+y^9+y^8+y^4+y^2+1)*(y^2+y+1)^3*(y^11+y^2+1)*x^13+(y^18+y^17+y^16+y^15+y^14+y^13+y^11+y^10+y^7+y^6+y^5+y^4+1)*(y^8+y^7+y^6+y^4+y^2+y+1)*(y^2+y+1)^2*x^12+(y^2+y+1)^3*(y^20+y^16+y^15+y^13+y^12+y^8+y^4+y^3+y^2+y+1)*y^4*x^11+(y^2+y+1)^4*(y^3+y^2+1)^2*(y^8+y^4+y^3+y+1)^2*x^10+(y^2+y+1)^5*(y^19+y^18+y^17+y^16+y^15+y^11+y^10+y^8+y^6+y^5+y^2+y+1)*y*x^8+(y^2+y+1)*(y^3+y+1)*(y+1)^10*x^29+(y^4+y^3+1)^2*(y^2+y+1)^6*(y^5+y^3+y^2+y+1)^2*x^6+(y^2+y+1)^6*(y^17+y^16+y^14+y^12+y^11+y^10+y^8+y^7+y^6+y^5+y^3+y+1)*y*x^5+(y^2+y+1)^7*(y^10+y^8+y^6+y^5+1)*(y^6+y^5+y^2+y+1)*x^4+(y^2+y+1)^4*(y^22+y^21+y^19+y^18+y^16+y^14+y^13+y^12+y^11+y^9+y^8+y^6+y^4+y^2+1)*x^9+(y^9+y^8+y^7+y^6+y^4+y^3+1)^2*(y+1)^2*x^18+(y+1)^10*y*x^28+(y^2+y+1)^14*y^2*x^2+(y^2+y+1)*(y^28+y^24+y^21+y^20+y^17+y^16+y^14+y^10+y^8+y^7+y^4+y^3+1)*x^15+(y^2+y+1)^5*(y^20+y^16+y^15+y^14+y^12+y^10+y^8+y^3+1)*x^7+(y+1)^20*y^2*x^22+(y^8+y^7+y^5+y^4+y^3+y^2+1)*(y^2+y+1)^7*(y^4+y^3+y^2+y+1)*(y^4+y+1)*x^3+(y^9+y^8+y^7+y^6+y^4+y^3+1)^2*(y^2+y+1)*(y^3+y+1)*(y+1)^2*x^17+(y^2+y+1)*(y^9+y^5+1)*(y^7+y^5+y^4+y^3+1)*(y+1)^4*x^23+(y^2+y+1)*(y^5+y^4+y^3+y+1)*(y+1)^7*y*x^25+(y^7+y^3+1)*(y^14+y^12+y^10+y^9+y^7+y^6+y^3+y^2+1)*(y+1)*(y^6+y^4+y^2+y+1)*y*x^16 ;

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

function FunctionDegrees(divisor)
    constantField := ConstantField(FunctionField(divisor));
    space,map := RiemannRochSpace(divisor);
    return [Degree(map(i)) : i in space | map(i) notin constantField];
end function;
function Run(s)
   return {Sort(FunctionDegrees(D)) : D in s};                  // for the actual proof
   // return Dimension(RiemannRochSpace(SetToSequence(s)[1]));  // for testing the file
end function;

print("Part I of the computation for N = 57");

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

// [ 18, 9, 0, 0, 0, 2, 36, 18, 4, 216 ]

// There is one orbit of 18 rational places, and one orbit of 9 degree-2 places.
// The genus is 85.

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

// Degree >1 places: [10]
// Rat-places: [] (only multiplicities >1 are listed).
// Do six deg10 places at a time (there are 216 in total)
plc10 := Places(AFF,10);
Run({ cuspsum + &+plc10[ 6*i-5 .. 6*i ]  : i in [1 .. 216/6] });
// Estimate for this: about 15000 seconds.

// Degree >1 places: [9]
// Rat-places: [2] (only multiplicities >1 are listed).
// Do all 4 places of degree 9 at the same time:
Run({ cuspsum + plc1[1] + &+Places(AFF,9) });

// Degree >1 places: [8] OR [8,2]
// Rat-places:    [2, 2] OR [] (only multiplicities >1 are listed).
plc8 := Places(AFF,8);
// Lets run three of them at a time.
Run({ 2*cuspsum + plc1[1] + plc2[1] + &+plc8[ 3*i-2 .. 3*i ] : i in [1 .. 18/3] });

// Degree >1 places: [7] OR [7,2]
// Rat-places: [2, 2, 2] OR [2] (only multiplicities >1 are listed).
// Lets run four degree-7 places at the same time:
plc7 := Places(AFF,7);
Run({ 2*cuspsum + 2*plc1[1] + plc2[1] + &+plc7[ 4*i-3 .. 4*i ] : i in [1 .. 36/4] });

// Degree >1 places: [6,2] OR [6,2,2]
// Rat-places:       [2,2] OR [] (only multiplicities >1 are listed).
// Lets run both degree-6 places at the same time:
Run({ 2*cuspsum + plc1[1] + cuspsum2 + plc2[1] + &+Places(AFF,6) });

// Degree >1 places: [6]
// Rat-places:       [2,2,2,2] (only multiplicities >1 are listed).
Run({ 3*cuspsum + 2*plc1[1] + &+Places(AFF,6) });

print("Now all cases that have a place of degree > 2 have been covered.");

print("The remaining cases are divided according to how many degree-2");
print("places they have. This can be: 0, 1, 2, 3, 4, or 5  degree-2");
print("places, counting with multiplicity. The remainder of this file");
print("will cover 0 completely (i.e. only degree-1 places) and cover");
print("a subset of the rest.");

// First consider: Only rational-places.
// We consider all partitions of 16:  P16 = {[16], [1,15], [2,14], [1,1,14], ...}
// Since f^(-1)({infinity} contains at least 6 distinct rational points, discard those
// with fewer than 6 entries.  Next, in each partition, remove 1's 2's and 3's.  What is left is:
// L16 = {[4], [5], [6], [7], [8], [9], [10], [11], [4, 4], [4, 5], [4, 6], [4, 7], [4, 8],
//        [5, 5], [5, 6], [5, 7], [6, 6], [4, 4, 4], [4, 4, 5]}
Run({
	4*cuspsum + plc1[1],   // Has dim=2, covers: {[4, 4], [4, 5],  [4, 4, 4], [4, 4, 5]}
	3*cuspsum + 8*plc1[1], //            covers: {[4], [5], [6], [7], [8], [9], [10], [11]}
	3*cuspsum + 3*plc1[1] + 3*( &+plc1[2..6]   ), //
	3*cuspsum + 3*plc1[1] + 3*( &+plc1[7..10]  ), // These 4 cover:  {[4, 6], [4, 7], [4, 8],
	3*cuspsum + 3*plc1[1] + 3*( &+plc1[11..14] ), //                  [5, 5], [5, 6], [5, 7],
	3*cuspsum + 3*plc1[1] + 3*( &+plc1[15..18] )  //                  [6, 6]}
});

print("Of the remaining cases, we now treat cases for which:");
print("max(multiplicity of deg1 places) + max(multiplicity of deg2 places) <= 4");

Run({	3*cuspsum - 4*plc1[1] + 1*cuspsum2,
	2*cuspsum - 3*plc1[1] + 2*cuspsum2,
	1*cuspsum - 2*plc1[1] + 3*cuspsum2	});

print("Not all cases are covered yet, continued in X1_57_Part_II.m");
