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

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

// Say f is a function of degree <= 8, since f^(-1)({0,1,2,infinity}) contains
// all 22 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 2.

//[ 0, 0, 0, -2, -8, -7, 0, -3, -3, 0, 2, 2, -1, -1, 1, 2, -1, -2, 0, 3, 8, 8, 2, 0]
// C0,C1,C2, C3, C4, C5,

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

// Contains one member of each diamond-orbit:
orbs := [orb1, orb2];

// This covers all cases:
{FunctionDegrees( cuspsum + orb1 + orb2 + P ) : P in plc1};

