// Is the degree 9 point on X_1(28) P^1 isolated?

E := EllipticCurve([1,-1,1,137,2643]);

// This is an elliptic curve with j-invariant 3^3 * 13/2^2.

R<x> := PolynomialRing(Rationals());

// Let's build the extension over which the curve has a point of order 28.
div2 := DivisionPolynomial(E,2);
div7 := DivisionPolynomial(E,7);

f1 := Factorization(div7)[1][1];
K1 := NumberField(f1);

// f1 has degree 3. We build the degree 3 number field.

SS := PolynomialRing(K1);
K0 := AbsoluteField(NumberField(SS!div2));

// We adjoin an x-coordinate of a 2-torsion point to K1.

OK := LLL(MaximalOrder(K0));
poly0 := MinimalPolynomial(OK.2);

// We choose a simple polynomial defining K.

poly := x^9 + x^8 - 2*x^7 + x^6 + 5*x^5 - x^4 - 5*x^3 - 2*x^2 + 2*x + 1;
IsIsomorphic(NumberField(poly0),NumberField(poly));

K<z> := NumberField(poly);
EK := BaseChange(E,K);

T, Tmap := TorsionSubgroup(EK);

// We see that over K, the torsion subgroup of E is Z/28Z.

Px := -182*z^8 + 65*z^7 + 260*z^6 - 546*z^5 - 156*z^4 + 364*z^3 + 364*z^2 - 156*z - 120;
Py := -5148*z^8 + 2080*z^7 + 7306*z^6 - 15444*z^5 - 3978*z^4 + 10634*z^3 + 
    10634*z^2 - 4654*z - 3574;

// These are the x and y coordinates of one point of order 28. There are
// 11 other points of order 28 in E(K).

// Drew Sutherland's model of X_1(28) has the form
// E_u : y^2 + xy + uy = x^3 + ux^2. This curve has j-invariant
// (-256u^6 + 768u^5 - 816u^4 + 352u^3 - 51u^2 + 3u - 1/16)/(u^5 - (1/16)u^4).
// We find the values of u that make this expression equal 351/4.

pol := (-256*x^6 + 768*x^5 - 816*x^4 + 352*x^3 - 51*x^2 + 3*x - 1/16)
- (351/4)*(x^5 - (1/16)*x^4);

uvals := Roots(PolynomialRing(K)!pol);
// There are two values for u in K that make j(E_u) = 351/4.

A<x,y> := PolynomialRing(K,2);
X := x^2*(x + 1)^3*y^6 
        + x^2*(x + 1)^2*(x^4 + 11)*y^5
        + 3*x^2*(x + 1)*(3*x^4 - 2*x^2 + 15)*y^4
        - (5*x^8 - 28*x^6 + 42*x^4 - 84*x^2 + 1)*y^3
        - (x - 1)*(25*x^6 - 23*x^4 + 67*x^2 - 5)*y^2
        + 6*(x - 1)^5*(x + 1)^3*y
        - (x - 1)^6*(x + 1)^3;

// This curve X is Drew Sutherland's model for X_1(28).

C<x,y,zz> := ProjectiveClosure(Curve(AffineSpace(K,2),X));
P1 := ProjectiveSpace(Rationals(),1);
q := (x*y + y*zz + 2*zz^2)/(2*x*zz);
t := (-x*y - y*zz - 2*zz^2)/(2*zz^2);

// According to https://math.mit.edu/~drew/X1/X1opt28.txt,
// the elliptic curve corresponding to a rational point on X_1(28)
// is [1,(q^2-1)*(t^2-1)/16,(q^2-1)*(t^2-1)/16,0,0]. The point is
// [(q+1)*(t^2-1)/8,(q+1)^2*(t-1)^2*(t+1)/32].

elt := (q^2-1)*(t^2-1)/16;

des1 := DefiningPolynomial(C);
des2 := Numerator(elt)-uvals[1][1]*Denominator(elt);

X1 := Scheme(ProjectiveSpace(K,2),[des1,des2]);

des2 := Numerator(elt)-uvals[2][1]*Denominator(elt);

X2 := Scheme(ProjectiveSpace(K,2),[des1,des2]);

// We find the points on X so that u = (q^2-1)*(t^2-1)/16 equals
// the two u-values found above.

// The RationalPoints(X1); command takes about 5 minutes. It only returns
// (1 : 0 : 0) and (0 : 1 : 0). These points are both singular, but
// all the points above them on the non-singular model are cusps.

// The RationalPoints(X2); command also takes a while. It returns the
// two points (1 : 0 : 0) and (0 : 1 : 0) found above. It also returns
// the points listed below.

pointlist := [ C(K)![-3*z^8 + 4*z^6 - 8*z^5 - 5*z^4 + 4*z^3 + 5*z^2 - 1,
3*z^8 - z^7 - 5*z^6 + 8*z^5 + 5*z^4 - 9*z^3 - 8*z^2 + 3*z + 3, 1],
C(K)![-3*z^8 + 8*z^6 - 8*z^5 - 9*z^4 + 14*z^3 + 11*z^2 - 2*z - 9,
z^8 + 2*z^7 + z^6 + 5*z^4 + 8*z^3 - 2*z^2 - 4*z - 3],C(K)![-z^8 - 2*z^5 - z^4 - 2*z^3 + z^2 + 1,-3*z^8 - 3*z^7 + 4*z^6 - 4*z^5 - 13*z^4 - z^3 + 9*z^2 + 5*z - 2],C(K)![z^8 + 2*z^5 + z^4 + 2*z^3 - z^2 - 1,9*z^8 + 5*z^7 - 20*z^6 + 18*z^5 + 37*z^4 - 25*z^3 - 33*z^2 - 3*z + 18],C(K)![3*z^8 - 8*z^6 + 8*z^5 + 9*z^4 - 14*z^3 - 11*z^2 + 2*z + 9,-3*z^8 + 2*z^7 + 3*z^6 - 10*z^5 + z^4 + 4*z^3 + 4*z^2 - 4*z - 1],C(K)![3*z^8 - 4*z^6 + 8*z^5 + 5*z^4 - 4*z^3 - 5*z^2 + 1,z^8 + z^7 - z^6 + 2*z^5 + 3*z^4 + z^3 - 3*z - 1]];

// There are six such points, and that's what we should expect. The points
// <E,P> and <E,-P> are identified, and since E has 12 points of order 28 on
// it, there should be 6 points on X_1(28) lying above E.

// The Riemann-Roch calculation is too difficult to do over K.
// Instead, we'll do it modulo p for some prime p.

p1 := MinimalPolynomial(pointlist[6][1]);
p2 := MinimalPolynomial(pointlist[6][2]);

p := 11;
F := GF(p);
R<x,y,zz> := PolynomialRing(F,3);

// Now, let's rebuild X_1(28) in characteristic p.

newp1 := &+[ Coefficient(p1,i)*x^i*zz^(9-i) : i in [0..9]];
newp2 := &+[ Coefficient(p2,i)*y^i*zz^(9-i) : i in [0..9]];

cpoly := x^2*(x + zz)^3*y^6*zz^2
        + x^2*(x + zz)^2*(x^4 + 11*zz^4)*y^5
        + 3*x^2*(x + zz)*(3*x^4 - 2*x^2*zz^2 + 15*zz^4)*y^4*zz^2
        - (5*x^8 - 28*x^6*zz^2 + 42*x^4*zz^4 - 84*x^2*zz^6 + zz^8)*y^3*zz^2
        - (x - zz)*(25*x^6 - 23*x^4*zz^2 + 67*x^2*zz^4 - 5*zz^6)*y^2*zz^4
        + 6*(x - zz)^5*(x + zz)^3*y*zz^4
        - (x - zz)^6*(x + zz)^3*zz^4;
C := Curve(ProjectiveSpace(F,2),cpoly);
XX := Scheme(AmbientSpace(C),[newp1,newp2,cpoly]);

// The subscheme XX consists of the points on C whose x and y coordinates
// satisfy the equations the degree 9 point satisfies. We use this to make
// a degree 9 divisor on C.

t0 := Cputime();
printf "Making divisor from point 1.\n";
D := Divisor(C,XX);
printf "Time taken was %o.\n",Cputime(t0);
printf "Degree of divisor is %o.\n",Degree(D);

printf "Computing Riemann-Roch space.\n";
V := RiemannRochSpace(D);
printf "Riemann-Roch space has dimension %o.\n",Dimension(V);

// If the Riemann-Roch space has dimension 1, that means that the only functions
// on C/F_p that have poles only at D are the constant functions.
// If f : C/K -> P^1/K is a non-constant function with poles only at D,
// then there must be alpha and beta in Q so that alpha*f + beta is non-constant
// mod p and the reduction of alpha*f + beta is non-constant. This is a
// contradiction, and so there is no non-constant f : C/K -> P^1/K with
// poles only at D.