A<x,y> := PolynomialRing(GF(3),2);
P,homogenize := Homogenization(A);
print("f := The polynomial from math.mit.edu/~drew/FFFc26.txt");
f := x^6*y^2 - 5*x^5*y^3 + 2*x^5*y^2 + x^4*y^5 + 5*x^4*y^4 - 2*x^4*y^2 - x^4*y -    2*x^3*y^6 - 4*x^3*y^4 + 9*x^3*y^2 - 4*x^3*y + x^2*y^7 + 5*x^2*y^4 - 11*x^2*y^3 + 4*x^2*y^2 + x^2*y - x*y^7 + 2*x*y^6 - 4*x*y^5 + 6*x*y^4 - 3*x*y^3 + x*y^2 - 2*x*y + x + y^6 - 3*y^5 + 3*y^4 - y^3;
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, 2 is:");
[#Places(AFF,i) : i in [1..2]];
print("There are no places of degree 2.");
print("");
print("12/4 = 3. So to show gon_{F_3} > 5 it suffices to look at divisors");
print("of the form cuspsum + D with deg D = 2.");
print("Both type(D)=(2) and type(D)=(1,1) follow from the following calculation:");
Min(&cat[FunctionDegrees(cuspsum + q + r) : q,r in div1]);
//
// Note: div1 is divided into two orbits under diamond operators.  Instead of
// looping q over div1, it would be sufficient to let q loop over a set of two
// elements, one from each orbit.
// That would be 6 times faster, but it would take a little bit more programming.
