clear;
q:=RandomPrime(50);
Fq:=GF(q);
Poly<x>:=PolynomialRing(Fq);
f4:=Random(Fq); f3:=Random(Fq); f2:=Random(Fq); f1:=Random(Fq); f0:=Random(Fq);
f:=x^5+f4*x^4+f3*x^3+f2*x^2+f1*x+f0;
C2:=HyperellipticCurve(f);
F<x,y>:=FunctionField(C2);
O:=PointsAtInfinity(C2)[1];

P1:=Random(C2); P2:=Random(C2); P3:=Random(C2); P4:=Random(C2);
D:=Divisor(P1)+Divisor(P2)+Divisor(P3)+Divisor(P4)-4*Divisor(O);

l:=F!Interpolation([P1[1],P2[1],P3[1],P4[1]],[P1[2],P2[2],P3[2],P4[2]]);
l:=F!(y-l);
Support(Divisor(l));

