A<x,y> := PolynomialRing(GF(2),2);
P,homogenize := Homogenization(A);
print("f := The polynomial from math.mit.edu/~drew/FFFc33.txt");
f := y^12 + (3*x^3 - 8*x^2 - 7*x - 7)*y^11 + (3*x^6 - 15*x^5 + 7*x^4 + 23*x^3 + 75*x^2 + 43*x + 21)*y^10 
          + (x^9 - 10*x^8 + 14*x^7 + 14*x^6 + 38*x^5 - 152*x^4 - 228*x^3 - 268*x^2 - 111*x - 35)*y^9
          - (2*x^11 - 8*x^10 - 6*x^9 - 12*x^8 + 110*x^7 + 2*x^6 - 189*x^5 - 696*x^4 - 677*x^3 - 496*x^2 - 155*x - 35)*y^8 
          + (x^13 - x^12 - 5*x^11 - 33*x^10 + 18*x^9 + 81*x^8 + 202*x^7 - 359*x^6 - 915*x^5 - 1428*x^4 - 1000*x^3 - 525*x^2 - 125*x - 21)*y^7
          + (4*x^12 + 21*x^11 + 24*x^10 - 134*x^9 - 174*x^8 + 39*x^7 + 1063*x^6 + 1550*x^5 + 1565*x^4 + 817*x^3 + 323*x^2 + 57*x + 7)*y^6
          - (26*x^11 - 38*x^10 - 215*x^9 - 36*x^8 + 484*x^7 + 1335*x^6 + 1323*x^5 + 957*x^4 + 368*x^3 + 110*x^2 + 13*x + 1)*y^5
          + (21*x^11 - 74*x^10 - 133*x^9 + 180*x^8 + 570*x^7 + 866*x^6 + 600*x^5 + 317*x^4 + 83*x^3 + 18*x^2 + x)*y^4 
          - (7*x^11 - 52*x^10 - 15*x^9 + 191*x^8 + 297*x^7 + 296*x^6 + 136*x^5 + 51*x^4 + 7*x^3 + x^2)*y^3
          + (x^11 - 16*x^10 + 19*x^9 + 80*x^8 + 73*x^7 + 49*x^6 + 12*x^5 + 3*x^4)*y^2 + (2*x^10 - 8*x^9 - 15*x^8 - 7*x^7 - 3*x^6)*y + x^9 + x^8;
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-4 is:");
[#Places(AFF,i) : i in [1..4]];     
print("There are 10 places of degree 1 and no places of degree 3 and 4");
print("");
print("ceil(10/3)=4 so to show gon_{F_2} > 9 it suffices to look at divisors");
print("of the form cuspsum + D with deg D = 5.");
print("");
print("The diamond operators act transitively on the places of degree 1, so");
print("instead of looping all variables over div1, we select one to be div1[1]");
print("saving a factor 10 in CPU time.");
print("");

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

print("type(D) = (3,1,1), (2,1,1,1):");
FunctionDegrees(2*cuspsum + 2*div1[1]);

print("type(D) = ((2,1),(1,3)), ((2,1),(1,2),(1,1)) and ((2,1),3(1,1)):");
Max([Dimension(RiemannRochSpace(cuspsum + p + q + r + div1[1])) : p in Places(AFF,2), q,r in div1]); 

print("type(D) = ((2,2),(1,1)) and ((2,1),(2,1),(1,1)):");
Max([Dimension(RiemannRochSpace(cuspsum + p + q + div1[1])) : p,q in Places(AFF,2)]); 

print("type(D) = ((5,1)):");
Max([Dimension(RiemannRochSpace(cuspsum+p)) : p in Places(AFF,5)]); 
