A<x,y> := PolynomialRing(GF(5),2);
P,homogenize := Homogenization(A);
// print("f := The polynomial from math.mit.edu/~drew/FFFc48.txt reduced mod 5");
// f := (y^2+4*y+1)*y^6*x^22+4*y^6*(y^3+4*y^2+y+2)*x^21+2*(y+3)*y^5*x^20+2*(y^3+4*y^2+3*y+4)*y^5*x^19+4*(y^3+3*y^2+4*y+1)*y^4*(y+2)*(y+4)*x^18+3*y^4*(y^2+y+2)*(y^7+3*y^6+y^5+3*y^4+3*y^3+y^2+y+1)*x^17+2*y^3*(y^4+y^3+4)*(y^5+3*y^4+y^3+3*y^2+3)*(y^2+2)*x^16+y^3*(y^8+3*y^7+4*y^6+4*y^5+y+1)*x^15+2*y^2*(y^9+4*y^8+4*y^7+2*y^6+y^5+y^4+3*y^3+3*y^2+4*y+1)*(y^2+y+2)*x^14+(y^3+4*y^2+4)*(y+3)^2*y^2*(y^3+3*y+2)*(y+4)*(y^3+4*y^2+y+1)*x^13+(y+3)^2*y*(y^5+3*y^3+2*y^2+y+1)*(y+1)*(y^9+2*y^8+4*y^7+4*y^6+2*y^5+y^4+3*y^3+4*y^2+4*y+1)*x^12+4*y*(y^15+2*y^14+4*y^12+4*y^11+y^10+y^9+y^7+4*y^6+3*y^5+4*y^4+3*y^3+y^2+4*y+4)*(y+2)*(y^2+2)*x^11+2*(y^5+3*y^2+y+3)*(y^2+3*y+4)*(y^2+4*y+1)*(y^6+y^5+4*y^4+3*y^3+y^2+1)*y*x^10+y*(y^7+2*y^6+3*y^5+y^4+y^3+3*y^2+4*y+2)*(y^10+y^8+y^7+2*y^6+2*y^5+2*y^4+y^3+2)*x^9+(y^5+4*y^4+4*y^3+4*y^2+3*y+3)*y*(y^10+2*y^9+3*y^8+y^7+y^6+y^5+4*y^4+2*y^3+3*y^2+3*y+1)*(y+2)*(y+4)*x^8+2*(y^7+4*y^6+y^5+3*y^4+3*y^3+y^2+y+2)*(y^3+4*y^2+3*y+4)*y*(y+1)*(y^3+4*y^2+4*y+4)*(y+4)^2*x^7+(y^5+3*y^4+y^3+2*y+2)*y*(y^10+y^8+4*y^7+2*y^6+2*y^5+3*y^4+y^3+4*y^2+3)*(y+4)^3*x^6+(y^6+4*y^5+y^4+2*y^3+4*y^2+4*y+3)*y*(y^3+y+4)*(y+4)^4*(y^4+4*y^3+3*y^2+y+2)*x^5+y*(y+2)*(y^10+2*y^9+2*y^8+2*y^7+y^6+y^5+3*y^4+3*y^2+3*y+4)*(y+4)^5*x^4+(y^8+y^6+4*y^5+3*y^4+4*y^3+3*y^2+y+4)*y^2*(y+4)^6*x^3+(y+4)^7*(y^8+y^7+2*y^6+2*y^5+y^3+3*y^2+1)*x^2+4*(y^2+y+1)*y^2*(y+2)*(y+4)^8*x+y^2*(y+4)^9 ;

print("f := An optimized equation for X1(48) reduced mod 5");
f := (y+3)*(y+4)*(y+1)^5*x^16+3*(y+3)*(y+4)*(y+1)^6*x^15+4*(y+1)^6*(y^3+3*y^2+2*y+3)*x^14+4*(y^2+4*y+1)*(y+1)^5*(y^3+y^2+1)*x^13+4*(y^4+y^3+2*y^2+3*y+4)*(y^6+2*y^4+y^3+3*y^2+2)*(y+1)^2*x^12+2*(y+3)*(y+1)^3*(y^9+2*y^8+4*y^7+4*y^6+3*y^5+2*y^4+4*y^3+y^2+2*y+3)*x^11+2*(y^12+4*y^9+y^8+y^7+3*y^4+3*y+1)*(y+4)*(y+1)*x^10+4*(y+3)^2*(y+1)^2*(y^11+4*y^10+4*y^9+3*y^8+y^7+4*y^6+2*y^5+y^4+2*y^3+2*y^2+3*y+2)*x^9+2*(y+3)*(y^6+4*y^5+y^4+4*y+2)*(y^3+2*y+4)*(y^5+3*y^4+4*y^3+y^2+y+3)*x^8+4*(y+3)*(y^3+2*y^2+1)*(y+1)*(y^12+y^11+4*y^10+3*y^9+2*y^8+3*y^7+4*y^6+2*y^5+y^4+4*y^3+2*y^2+y+1)*x^7+2*(y^3+y+1)*(y^3+2*y^2+1)*(y+1)^2*(y^10+y^9+2*y^8+4*y^7+y^6+y^5+y^4+4*y+4)*x^6+2*(y^11+3*y^10+2*y^9+3*y^8+2*y^7+4*y^6+2*y^5+3*y^4+4*y^3+4)*(y+1)^3*(y^5+3*y^3+y^2+1)*x^5+4*(y^10+y^9+4*y^8+2*y^6+2*y^5+3*y^4+y^3+4*y^2+y+2)*(y+1)^4*(y^2+y+1)^2*(y+2)^2*x^4+3*(y^6+3*y^5+2*y^4+2*y^3+y^2+4*y+3)*(y+1)^6*(y^8+y^7+y^6+y^5+2*y^4+y^2+1)*x^3+(y+1)^7*(y^14+4*y^13+y^12+2*y^11+4*y^10+y^9+2*y^8+2*y^6+2*y^5+y^3+3*y^2+4*y+2)*x^2+3*(y+1)^13*(y^2+2*y+4)*(y^2+2*y+3)*x+(y+1)^14*(y^2+2*y+4)*(y^2+2*y+3);

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

function FunctionDegrees(divisor)
    constantField := ConstantField(FunctionField(divisor));
    space,map := RiemannRochSpace(divisor);
    return [Degree(map(i)) : i in space | map(i) notin constantField];
end function;

print("The number of places of degrees 1..5 is:");
[#Places(AFF,i) : i in [1..5]];
// [ 16, 12, 16, 68, 656 ]

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

// Contains one member of each diamond-orbit:
orbs := plc1;


// The cases below are semi-automatically generated:

// Degree >1 places: [2, 3]
// Rat-places: [] (only multiplicities >1 are listed).
s := {(cuspsum + P1+P2) : P1 in Places(AFF,2), P2 in Places(AFF,3)};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;

// Degree >1 places: [5]
// Rat-places: [] (only multiplicities >1 are listed).
s := {(cuspsum + P1) : P1 in Places(AFF,5)};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;

// Degree >1 places: [4]
// Rat-places: [2] (only multiplicities >1 are listed).
s := {(cuspsum + P1 + 1*Q1) : P1 in Places(AFF,4), Q1 in orbs};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;

// Degree >1 places: [2, 2]
// Rat-places: [2] (only multiplicities >1 are listed).
s := {(cuspsum + P1+P2 + 1*Q1) : P1 in Places(AFF,2), P2 in Places(AFF,2), Q1 in orbs};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;

// Degree >1 places: [3]
// Rat-places: [2, 2] (only multiplicities >1 are listed).
s := {(cuspsum + P1 + 1*Q1+1*Q2) : P1 in Places(AFF,3), Q1 in plc1, Q2 in orbs};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;

// Degree >1 places: [2]
// Rat-places: [2, 2, 2] (only multiplicities >1 are listed).
s := {(cuspsum + P1 + 1*Q1+1*Q2+1*Q3) : P1 in Places(AFF,2), Q1 in plc1, Q2 in plc1, Q3 in orbs};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;

// Degree >1 places: []
// Rat-places: [2, 2, 2, 3] (only multiplicities >1 are listed).
s := {(cuspsum + 1*Q1+1*Q2+1*Q3+2*Q4) : Q1 in plc1, Q2 in plc1, Q3 in plc1, Q4 in orbs};
print "The number of RR-computations to be done for this case is:", #s;
s := &cat[FunctionDegrees(D) : D in s];
if not s eq [] then s := Min(s); end if;
s;
