A<y,x> := PolynomialRing(GF(3),2);
P,homogenize := Homogenization(A);
print("f := An optimized equation for X1(44) reduced mod 3");
f := (y+1)^7*x^15+(y+1)^6*(y^3+2*y^2+2*y+2)*(y+2)*x^12+(y+1)^7*y*(y^2+2*y+2)*(y+2)*x^11+2*(y+1)^4*y*(y^2+2*y+2)*(y+2)*(y^4+2*y^3+y+1)*x^10+2*(y^2+1)*(y+1)^4*(y^5+2*y^4+1)*(y+2)*x^9+2*(y^2+1)*(y^2+y+2)*(y+1)^2*y*(y^2+2*y+2)*(y+2)*x^8+2*(y^5+y^4+y^2+1)*(y+1)^2*y*(y^2+2*y+2)*(y+2)*x^7+2*(y^10+2*y^8+y^2+1)*y*x^6+(y^2+1)*(y^2+y+2)*(y+1)^2*y*(y^2+2*y+2)*(y+2)^4*x^5+(y^5+y^4+y^2+1)*(y+1)^2*y*(y^2+2*y+2)*(y+2)^4*x^4+2*(y^9+2*y^8+y^6+y^4+2*y^3+2*y+2)*y^3*(y+2)^3*x^3+(y+1)^2*y^6*(y^2+2*y+2)*(y+2)^3*x^2+(y+1)*y^6*(y^2+2*y+2)*(y+2)^5*x+2*y^8*(y+2)^7 ;

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..4 is:");
[#Places(AFF,i) : i in [1..4]];
// [ 15, 5, 0, 0 ]

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

//[ 0, 0, 0, -3, 0, 0, -3, -4, 0, 0, 0, 0, 0, -1, -4, 2, 0, -2, 7, 4, 0, 4, 0 ]
// C0,C1,                                            C15       C18

pl,val := Support(Divisor(AFF.1));
for i in [1..#val] do
      if val[i] eq 2 then
                orb1 := pl[i];  // Cusp number 15
      end if;
      if val[i] eq 7 then
		orb2 := pl[i];  // Cusp number 18
      end if;
end for;

// Contains one member of each diamond-orbit:
orbs := [orb1, orb2];
// orb1 is in an orbit of length 10 and orb2 is in an orbit of length 5


// Covers the case of 4 rational points, each with multiplicity 2.
FunctionDegrees( 2*cuspsum );

// Degree >1 places: [2, 2]
// Rat-places: [] (only multiplicities >1 are listed).
s := {(cuspsum + P1+P2) : P1 in Places(AFF,2), P2 in Places(AFF,2)};
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] (only multiplicities >1 are listed).
s := {(cuspsum + P1 + 1*Q1+1*Q2) : P1 in Places(AFF,2), 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: []
// Rat-places: [2, 2, 3] (only multiplicities >1 are listed).
s := {(cuspsum + 1*Q1+1*Q2+2*Q3) : 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;
