A<x,y> := PolynomialRing(GF(3),2);
P,homogenize := Homogenization(A);
print("f := The polynomial from math.mit.edu/~drew/FFFc34.txt");
f := x^12*y^4 + x^11*y^6 - 10*x^11*y^5 + 2*x^11*y^4 - 5*x^10*y^7 + 33*x^10*y^6 - 6*x^10*y^4 - x^10*y^3 + 10*x^9*y^8 - 50*x^9*y^7 - 29*x^9*y^6 + 21*x^9*y^5 + 16*x^9*y^4 - 3*x^9*y^3 - 10*x^8*y^9 + 35*x^8*y^8 + 75*x^8*y^7 - 20*x^8*y^6 - 41*x^8*y^5 - 18*x^8*y^4 + 14*x^8*y^3 + 5*x^7*y^10 - 6*x^7*y^9 - 80*x^7*y^8 - 10*x^7*y^7 + 25*x^7*y^6 + 90*x^7*y^5 - 35*x^7*y^4 - 10*x^7*y^3 - x^6*y^11 - 5*x^6*y^10 + 38*x^6*y^9 + 30*x^6*y^8 + 25*x^6*y^7 - 120*x^6*y^6 + 10*x^6*y^5 + 6*x^6*y^4 + 39*x^6*y^3 - 15*x^6*y^2 + 2*x^5*y^11 - 5*x^5*y^10 - 19*x^5*y^9 - 34*x^5*y^8 + 45*x^5*y^7 + 40*x^5*y^6 + 33*x^5*y^5 - 92*x^5*y^4 + 23*x^5*y^3 + 6*x^5*y^2 - x^4*y^11 + 4*x^4*y^10 + 9*x^4*y^9 + 18*x^4*y^8 - 40*x^4*y^7 - 34*x^4*y^6 + 37*x^4*y^5 + 25*x^4*y^4 - 8*x^4*y^3 - 16*x^4*y^2 + 6*x^4*y + x^3*y^10 - 12*x^3*y^9 + 11*x^3*y^8 - 5*x^3*y^7 + 46*x^3*y^6 - 64*x^3*y^5 + 17*x^3*y^4 + 7*x^3*y^3 + x^3*y^2 - 2*x^3*y + 10*x^2*y^8 - 30*x^2*y^7 + 36*x^2*y^6 - 36*x^2*y^5 + 43*x^2*y^4 - 32*x^2*y^3 + 8*x^2*y^2 + 2*x^2*y - x^2 - 5*x*y^7 + 21*x*y^6 - 35*x*y^5 + 30*x*y^4 - 15*x*y^3 + 5*x*y^2 - x*y + y^6 - 5*y^5 + 10*y^4 - 10*y^3 + 5*y^2 - y; 
f := homogenize(f);
XF2 := Scheme(ProjectiveSpace(P),f);
FF := FunctionField(XF2);
AFF := AlgorithmicFunctionField(FF);
plc1 := Places(AFF,1);
divgrp := DivisorGroup(AFF);
div1 := [divgrp ! i : i in plc1];
cuspsum := &+ div1;
load "function_degrees.m";

print("The number of places of degrees 1-5:");
[#Places(AFF,i) : i in [1..5]];
print("There are 16 places of degree 1 and no places of degrees 2 and 3.");

// 16/4=4 so to show gon_{F_3} > 9 it suffices to look at divisors
// of the form cuspsum + D with deg D = 5

// Here we take one element from each of the two orbits under
// the diamond operators.  This way, the first two searches below
// save a factor 8 in CPU time  (they still take a while though).
//
pl,val := Support(Divisor(AFF.1));
for i in [1..#val] do
	if val[i] eq -3 then
		orb1 := pl[i];
	end if;
	if val[i] eq -11 then
		orb2 := pl[i];
	end if;
end for;

print("type(D) = (5), (4,1), (3,2), (3,1,1), (2,2,1) and (2,1,1,1):");
Min(&cat[FunctionDegrees(cuspsum + 2*p + q + r + s) : q,r,s in div1, p in [orb1,orb2]]); 

// type(D) = (1,1,1,1,1) not needed because a function in this set and not
// in one of the previous ones would have to have poles of order 2 at 5 distinct cusps
// and hence have degree at least 10.

print("type(D) = ((4,1),(1,1)):");
Min(&cat[FunctionDegrees(cuspsum + p + q) : p in Places(AFF,4), q in [orb1,orb2]]);

print("type(D) = ((5,1)):");
Min(&cat[FunctionDegrees(cuspsum + p) : p in Places(AFF,5)]);
