Q := RationalField();
F<b,c> := FunctionField(Q,2);
E := EllipticCurve([F | 1-c, -b, -b, 0, 0]);
EP := E![0,0];

//E is an elliptic curve in Tate normal form, with point EP=(0,0).

P<x,y> := PolynomialRing(Q,2);
j := jInvariant(E);
Nj := Numerator(F!j);
Dj := Denominator(F!j);
PJSeq := [P!Nj+12288000*P!Dj];
x1:=(13*EP)[1];
x2:=(14*EP)[1];
f27:=P!(Numerator(x1)*Denominator(x2)-Numerator(x2)*Denominator(x1));

//f27 is a function of b,c such that our elliptic curve EP has j-invariant -2^15*3*5^3 = -12,288,000 and the point (0,0) has order dividing 27.

F27:=Factorization(f27)[3][1];

//We take the third polynomial in the factorization to get order 27, rather than 3 or 9.

A<b,c> := PolynomialRing(Rationals(),2);
I:=ideal<A|PJSeq,F27>;
P:=PrimaryDecomposition(I);
I1:=ideal<A|P[1]>;
RadicalDecomposition(I1);

//We create and decompose ideals to find a Grobner basis.

B<T> := PolynomialRing(Rationals());
K<T> := NumberField(T^9-153*T^8+4788*T^7-22452*T^6-322146*T^5-75726*T^4-1116*T^3+540*T^2+9*T-1);
b:=((-176330769/55338945968)*442711567744*T^8+215916859239*T^7-6767659352433*T^6+32096518996263*T^5+452352147208959*T^4+78995182405613*T^3+2049024574405*T^2+148357480453*T-12095294539)/442711567744;
c:= T;
E1:= EllipticCurve([1-c,-b,-b,0,0]);

//K, b, and c are defined using the Grobner basis so that E1 is an elliptic curve defined over a field of degree 9.

TorsionSubgroup(E1);

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

r:= b/T;
s:= T/((b/T)-1);
xpoint:= (s-r)/(r*s-2*r+1);
ypoint:= (r*s-2*r+1)/(s^2-s-r+1);

//xpoint and ypoint are defined by the inverse maps in Drew Sutherland's "alternate" model of X_1(N) at <https://math.mit.edu/~drew/X1_altcurves.html>.

p1:= MinimalPolynomial(xpoint);
p2:= MinimalPolynomial(ypoint);
F:= RationalField();
A<x,y,z>:= PolynomialRing(K,3);
newp1:= &+[ Coefficient(p1,i)*x^i*z^(9-i) : i in [0..9]];
newp2:= &+[ Coefficient(p2,i)*y^i*z^(9-i) : i in [0..9]];
cpoly:= x^10*y^3 - 6*x^9*y^4 + 15*x^8*y^5 + 3*x^8*y^3*z^2 - 3*x^8*y^2*z^3 - 19*x^7*y^6 - 
    6*x^7*y^5*z - 3*x^7*y^4*z^2 + 5*x^7*y^3*z^3 + 3*x^7*y^2*z^4 + 12*x^6*y^7 + 18*x^6*y^6*z - 
    6*x^6*y^5*z^2 + 3*x^6*y^4*z^3 - 12*x^6*y^3*z^4 - 3*x^6*y^2*z^5 + 3*x^6*y*z^6 - 3*x^5*y^8 - 
    18*x^5*y^7*z + 6*x^5*y^6*z^2 - 12*x^5*y^5*z^3 + 27*x^5*y^4*z^4 - 3*x^5*y^3*z^5 - x^5*y^2*z^6 - 
    2*x^5*y*z^7 + 6*x^4*y^8*z + 3*x^4*y^7*z^2 + 10*x^4*y^6*z^3 - 24*x^4*y^5*z^4 + 3*x^4*y^4*z^5 - 
    2*x^4*y^3*z^6 + 5*x^4*y^2*z^7 + x^4*y*z^8 - x^4*z^9 - 3*x^3*y^8*z^2 - 3*x^3*y^7*z^3 + 15*x^3*y^5*z^5 - 
    6*x^3*y^4*z^6 - 3*x^3*y^3*z^7 - x^3*y^2*z^8 + x^3*y*z^9 + 6*x^2*y^7*z^4 - 12*x^2*y^6*z^5 + 
    10*x^2*y^5*z^6 - 13*x^2*y^4*z^7 + 15*x^2*y^3*z^8 - 7*x^2*y^2*z^9 + x^2*y*z^10 - 4*x*y^6*z^6 + 
    13*x*y^5*z^7 - 16*x*y^4*z^8 + 10*x*y^3*z^9 - 4*x*y^2*z^10 + x*y*z^11 + y^5*z^8 - 4*y^4*z^9 + 6*y^3*z^10 - 4*y^2*z^11
    + y*z^12;

//The curve cpoly is Sutherland's "alternate" model of X_1(27), homogenized.

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.

D:=Divisor(C,XX);

//This is the most time-consuming part of the code.

V:=RiemannRochSpace(D);
Dimension(V);