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};
end function;

print("This file covers all partitions not yet covered in Parts I, II, III.");

// Three, four or five deg-2 places (counting with multiplicity).
// The multiplicities > 1 of the deg1 places and all multiplicities of
// the deg2 places are indicated as:
//
//            [[deg1 mult's], [deg2 mult's]]
//
// So the deg=1 mult=1 places are omitted as usual, which is OK provided
// that every divisor contains at least 1*cuspsum.

partA := &+plc2[1..4];
partB := &+plc2[5..9];

Run({
	2*cuspsum + 3*plc1[1] + 3*partA,       //  [[5], [3]]   and   [[2, 4], [3]]
	2*cuspsum + 3*plc1[1] + 3*partB,       //

	2*cuspsum + 2*cuspsum2 + 2*plc1[1],    // [[2, 4], [1, 2]] and [[2, 4], [1, 1, 1]]

	3*cuspsum + 3*plc2[1] + 1*partA,       // [[3], [4]] and   [[3], [1, 3]] and
	3*cuspsum + 3*plc2[1] + 1*partB,       // [[3, 3], [3]]    [[3, 3], [1,2]]
	                                       // [[2, 2, 3], [3]] [[2, 2, 3], [1, 2]]
	                                       // [[2, 2, 2, 2], [3]]

	1*cuspsum + 2*cuspsum2 + 4*plc1[1],    //  [[3], [2, 2]] and [[3], [1, 1, 2]]
	                                       //  [[5], [1, 2]] and [[5], [1, 1, 1]]

	2*cuspsum + 1*cuspsum2 + 4*plc2[1]     //  [[2, 2], [4]] and [[2, 2], [1, 3]]
});	                                       //  [[],     [5]] and [[],     [1, 4]]


// Note: cases where
//	max(multiplicity of deg1 places) + max(multiplicity of deg2 places) <= 4
// were already handled in Part I.

print("This completes the proof that gon_Q(X1(57)) > 16");
