// Instead of working with the curves themselves, we work with their function fields.
// R is the function field of P^1.

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

// We construct the function field of E as a degree-3 extension of R.

P<y> := PolynomialRing(R);
G<beta> := FunctionField(-y^3 - 2*y - 1 + x^2);

// The inclusion of R into G corresponds to the map theta.

// The function field of C is the compositum of the function fields
// of E and F. The map F to P^1 factors as [V: U: W] to [V: W], followed by the
// involution [V: W] to [W: V]. Thus, pulling back x to the function field of F
// gives 1/v, so that 1/x and u obey the equation for F. We adjoin an element u
// to the function field of E obeying such a relation.

T<u> := PolynomialRing(G);
H<gamma> := FunctionField(u^2-((1/x)^3 + 1/x + 1));

// The inclusion of G into H gives the map phi_C. We verify that H has
// has the correct genus, and hence is the function field of C.

Genus(H);

/*
7
*/

// We now repeat the construction over F_5 to give the reduction of C at 5.

R<x> := FunctionField(FiniteField(5));
P<y> := PolynomialRing(R);
G<beta> := FunctionField(-y^3 - 2*y - 1 + x^2);
T<u> := PolynomialRing(G);
H<gamma> := FunctionField(u^2-((1/x)^3 + 1/x + 1));

// We verify that the reduction has the correct genus, so that 5 is a prime of good reduction.

Genus(H);

/*
7
*/

// We compute the numerator of the zeta function at 5.

zeta5 := LPolynomial(H);

// We factor the numerator of the zeta function.
A<t> := Parent(zeta5);
Factorization(zeta5);

/*
[
    <5*t^2 + t + 1, 1>,
    <5*t^2 + 3*t + 1, 1>,
    <3125*t^10 + 1250*t^9 + 125*t^8 + 50*t^7 + 95*t^6 + 70*t^5 + 19*t^4 + 2*t^3 + t^2 + 2*t + 1, 1>
]
*/

// As there is an irreducible factor of degree 10, the Jacobian of C must have a simple
// factor over Q of dimension 5. We show that this simple factor is absolutely simple.
// The first step is to show that the irreducible factor of degree 10 satisfies
// the conditions of [Sto08, Lemma 3].

Q5 := Factorization(zeta5)[3][1];

// In order to calculate the resultant, we create a polynomial ring in two variables X and T and
// map our polynomial to a polynomial in T.

R<T, X> := PolynomialRing(Rationals(), 2);
coerce_map := hom<A->R |T>;
Q5M := coerce_map(Q5);

// We evaluate at X and at TX, and take the resultant with respect to X.

resultant_5 := Resultant(Evaluate(Q5M, [X,0]), Evaluate(Q5M, [T*X,0]), X);

// We now have a polynomial in a single variable T. We do the substitution T = t.

coerce_back_map := hom<R->A | t,0>;
res_5 := coerce_back_map(resultant_5);

// We divide the polynomial by (1-t)^10, and check that the remainder is 0.

quotient, remainder := Quotrem(res_5, (1-t)^10);
remainder;

/*
0
*/

// Factoring the quotient, we verify that no monic irreducible factor has constant term 1.

Factorization(quotient);

/*
[
    <5, 26>,
    <625*t^10 - 250*t^9 + 775*t^8 - 280*t^7 + 321*t^6 + 2*t^5 + 321*t^4 - 280*t^3 + 775*t^2 - 250*t + 625, 1>,
    <9765625*t^40 + 46875000*t^39 + 130859375*t^38 + 283437500*t^37 + 525078125*t^36 + 864012500*t^35 + 1298696250*t^34 
        + 1827300000*t^33 + 2446234750*t^32 + 3146731000*t^31 + 3923545850*t^30 + 4768590900*t^29 + 5663552795*t^28 +
        6574297380*t^27 + 7459339969*t^26 + 8275444872*t^25 + 8979954607*t^24 + 9541962952*t^23 + 9946943628*t^22 +
        10190134728*t^21 + 10271470452*t^20 + 10190134728*t^19 + 9946943628*t^18 + 9541962952*t^17 + 8979954607*t^16 +
        8275444872*t^15 + 7459339969*t^14 + 6574297380*t^13 + 5663552795*t^12 + 4768590900*t^11 + 3923545850*t^10 +
        3146731000*t^9 + 2446234750*t^8 + 1827300000*t^7 + 1298696250*t^6 + 864012500*t^5 + 525078125*t^4 +
        283437500*t^3 + 130859375*t^2 + 46875000*t + 9765625, 2>
]
*/

// Having verified the Q_5 satisfies the conditions of [Sto08, Lemma 3], we repeat our process modulo 7. 

R<x> := FunctionField(FiniteField(7));
P<y> := PolynomialRing(R);
G<beta> := FunctionField(-y^3 - 2*y - 1 + x^2);
T<u> := PolynomialRing(G);
H<gamma> := FunctionField(u^2-((1/x)^3 + 1/x + 1));

// We verify that the reduction has the correct genus, so that 7 is a prime of good reduction.

Genus(H);

/*
7
*/

// We compute the numerator of the zeta function at 7.

zeta7 := LPolynomial(H);

// We factor the numerator of the zeta function at 7.

A<t> := Parent(zeta7);
Factorization(zeta7);

/*
[
    <7*t^2 - 3*t + 1, 2>,
    <16807*t^10 + 3773*t^8 + 882*t^7 + 609*t^6 + 168*t^5 + 87*t^4 + 18*t^3 + 11*t^2 + 1, 1>
]
*/

// We pick out the degree-10 factor.

Q7 := Factorization(zeta7)[2][1];

// In order to calculate the resultant, we create a polynomial ring in two variables X and T and
// map our polynomial to a polynomial in T.

R<T, X> := PolynomialRing(Rationals(), 2);
coerce_map := hom<A->R | T>;
Q7M := coerce_map(Q7);

// We evaluate at X and at TX, and take the resultant with respect to X.

resultant_7 := Resultant(Evaluate(Q7M, [X,0]), Evaluate(Q7M, [T*X,0]), X);

// We now have a polynomial in a single variable T. We do the substitution T = t.

coerce_back_map := hom<R->A | t,0>;
res_7 := coerce_back_map(resultant_7);

// We divide the polynomial by (1-t)^10, and check that the remainder is 0.

quotient, remainder := Quotrem(res_7, (1-t)^10);
remainder;

/*
0
*/

// Factoring the quotient, we verify that no monic irreducible factor has constant term 1.

Factorization(quotient);

/*
[
    <7, 26>,
    <2401*t^10 + 7546*t^9 + 14455*t^8 + 19656*t^7 + 22465*t^6 + 23230*t^5 + 22465*t^4 + 19656*t^3 + 14455*t^2 + 7546*t +
        2401, 1>,
    <282475249*t^40 + 968486568*t^39 + 1516142663*t^38 + 1528495808*t^37 + 1351198765*t^36 + 1526008372*t^35 +
        2232185690*t^34 + 3090480764*t^33 + 3603495574*t^32 + 3709099100*t^31 + 3685473554*t^30 + 3660163780*t^29 +
        3725720691*t^28 + 4119561488*t^27 + 4842750729*t^26 + 5485599496*t^25 + 5664747807*t^24 + 5454030008*t^23 +
        5362444852*t^22 + 5784206920*t^21 + 6131187524*t^20 + 5784206920*t^19 + 5362444852*t^18 + 5454030008*t^17 +
        5664747807*t^16 + 5485599496*t^15 + 4842750729*t^14 + 4119561488*t^13 + 3725720691*t^12 + 3660163780*t^11 +
        3685473554*t^10 + 3709099100*t^9 + 3603495574*t^8 + 3090480764*t^7 + 2232185690*t^6 + 1526008372*t^5 +
        1351198765*t^4 + 1528495808*t^3 + 1516142663*t^2 + 968486568*t + 282475249, 2>
]
*/

// We now check that the number fields generated by Q_5(t) and Q_7(t) are linearly disjoint by checking
// that Q_5 remains irreducible over the number field generated by a root of Q_7. This suffices, if alpha is a root
// of Q_5, then the powers of alpha form a basis for the extension defined by Q_5, and this basis remains 
// linearly independent over the number field generated by a root of Q_7.

F1 := NumberField(Q7);
Q5 := ChangeRing(Q5, F1);
R<t> := Parent(Q5);
Factorization(Q5);

/*
[
    <t^10 + 2/5*t^9 + 1/25*t^8 + 2/125*t^7 + 19/625*t^6 + 14/625*t^5 + 19/3125*t^4 + 2/3125*t^3 + 1/3125*t^2 + 2/3125*t 
        + 1/3125, 1>
]
/*