// Magma code accompanying the paper
// "Rational 6-cycles under iteration of quadratic polynomials"
//
// Michael Stoll, 2008-03-19
//
// To go through the computations, use 
//  iload "Xdyn06.magma";
//
///////////////
// Section 2 //
///////////////
//
// Set up the iteration
P2xc<x,c> := PolynomialRing(Integers(), 2);
f := x^2 + c;
f2 := Evaluate(f, [f,c]);
f3 := Evaluate(f2, [f,c]);
f4 := Evaluate(f3, [f,c]);
f5 := Evaluate(f4, [f,c]);
f6 := Evaluate(f5, [f,c]);
//
// Find equation of Y_1^dyn(6)
g1 := f - x;
g2 := ExactQuotient(f2 - x, g1);
g3 := ExactQuotient(f3 - x, g1);
g6 := ExactQuotient(f6 - x, g1*g2*g3);
//
// Now compute an equation of Y_0^dyn(6)
tr := x + f + f2 + f3 + f4 + f5;
P3<tt,xx,cc> := PolynomialRing(Integers(), 3);
e1 := Resultant(Evaluate(g6, [xx,cc]), tt-Evaluate(tr, [xx,cc]), xx);
flag, e2 := IsPower(e1, 6); assert flag;
h6 := -e2; h6;
Coefficients(h6, cc);
//
// Find a nicer model
P2uv<u,v> := PolynomialRing(Rationals(), 2);
h6a := -Numerator(Evaluate(h6, [2/u-1, 0, v/4-1/u^2+2/u-u/4]))/32;
P2uw<u,w> := PolynomialRing(Rationals(), 2);
h6b := Numerator(Evaluate(h6a, 
                          [u, ((2*u^2+3*u+6)*w-18)/(-u^2*w-u^2+3*u)]));
facth6b := Factorization(h6b); facth6b;
g := facth6b[2,1]; Coefficients(g, u);
//
// Find primes of bad reduction
P2Zuw<U,W> := PolynomialRing(Integers(), 2);
G := P2Zuw!g;
// affine patches
G1 := &+[Coefficient(G, U, i)*U^(3-i) : i in [0..3]];
G2 := &+[Coefficient(G, W, i)*W^(3-i) : i in [0..3]];
G3 := &+[Coefficient(G1, W, i)*W^(3-i) : i in [0..3]];
// singularities
I0 := ideal<P2Zuw | G, Derivative(G, U), Derivative(G, W)>;
I1 := ideal<P2Zuw | G1, Derivative(G1, U), Derivative(G1, W)>;
I2 := ideal<P2Zuw | G2, Derivative(G2, U), Derivative(G2, W)>;
I3 := ideal<P2Zuw | G3, Derivative(G3, U), Derivative(G3, W)>;
[Basis(EliminationIdeal(I, {}))[1] : I in [I0, I1, I2, I3]];
bad := &join[Set(PrimeDivisors(Integers()!B)) : B in $1]; bad;
bigp := Max(bad); // the large bad prime
//
// Find singularities at bigp
I0p := ChangeRing(I0, GF(bigp));
GBp := GroebnerBasis(I0p); GBp;
sptp := Variety(I0p)[1];
Gp := Generic(I0p)!G;
Gps := Evaluate(Gp, [Parent(Gp).1+sptp[1], Parent(Gp).2+sptp[2]]);
assert HomogeneousComponent(Gps, 0) eq 0;
assert HomogeneousComponent(Gps, 1) eq 0;
Factorization(HomogeneousComponent(Gps, 2));
// ==> simple double point with F_p-rational tangent directions
// check for regularity:
Valuation(Evaluate(G, [Integers()!sptp[1], Integers()!sptp[2]]), bigp);
// constant term not divisible by p^2 ==> regular point
//
// Find singularities at 2
Is2 := [ChangeRing(I, GF(2)) : I in [I0, I1, I2, I3]];
GBs := [*GroebnerBasis(I) : I in Is2*]; GBs;
// the first ideal gives all the singularities
sings := Variety(Is2[1]); sings;
//
// Analysis of these singularities, minimal regular model
PF2uw<uu,ww> := PolynomialRing(GF(2), 2);
//
// First singularity: <0,1> mod 2
// Shift to origin (mod 2)
Gsing1a := Evaluate(G, [U, W+1]);
[HomogeneousComponent(Gsing1a, d) : d in [0,1]];
// This shows that the point (u = 0, w = 1, 2 = 0) is non-regular,
// so we need to blow it up.
//
// A helper function:
function ReplaceP(pol, prime, rep)
  replc := func<c | ExactQuotient(c, prime^v)*rep^v
                     where v := Valuation(c, prime)>;
  repeat
    old := pol;
    pol := &+[replc(cs[i])*ms[i] : i in [1..#ms]]
            where cs := Coefficients(old) where ms := Monomials(old);
  until pol eq old;
  return pol;
end function;
//
// Define polynomial ring in three variables
P3Zpxy<x,y,p> := PolynomialRing(Integers(), 3);
P3F2pxy<xx,yy,pp> := PolynomialRing(GF(2), 3);
//
// Replace all 2's in the coefficients by p's and u by x, v by y
Gsing1b := ReplaceP(Evaluate(Gsing1a, [x,y]), 2, p);
// First chart: (x, y, p) <-- (x, x*y, x*p)
Gsing1c := ExactQuotient(Evaluate(Gsing1b, [x,x*y,x*p]), x^2);
// Look at special fiber: p = 0 or x = 0
P3F2pxy!Evaluate(Gsing1c, [x,y,0]);
P3F2pxy!Evaluate(Gsing1c, [0,y,p]);
// The first is the original curve, the second the line y+p+1 = 0, with
// multiplicity 2. The two intersect at x = p = 0, y = 1.
// We change coordinates so that the line is y = x = 0.
Gsing1d := Evaluate(Gsing1c, [x, y+p+1, p]);
Gsing1d := ReplaceP(Gsing1d, 2, x*p);
Coefficient(Coefficient(Gsing1d, x, 0), y, 0);
Coefficient(Coefficient(Gsing1d, x, 0), y, 1);
Coefficient(Coefficient(Gsing1d, x, 1), y, 0);
// All are zero ==> all the points on the line are non-regular.
// Before we blow up the line, we look at the other charts.
//
// Second chart: (x, y, p) <-- (y*x, y, y*p)
Gsing1e := ExactQuotient(Evaluate(Gsing1b, [y*x,y,y*p]), y^2);
P3F2pxy!Evaluate(Gsing1e, [x,y,0]);
P3F2pxy!Evaluate(Gsing1e, [x,0,p]);
// This is just another view of the same configuration as before.
//
// Third chart: (x, y, p) <-- (p*x, p*y, p)
Gsing1f := ExactQuotient(Evaluate(Gsing1b, [p*x,p*y,p]), p^2);
// Here the special fiber is p = 0 (since p = 2).
P3F2pxy!Evaluate(Gsing1f, [x,y,0]);
// This only shows the double line
//
// We blow up the line in the first chart.
// First sub-chart: (x, y, p) <-- (x, x*y, p)
Gsing1g := ExactQuotient(Evaluate(Gsing1d, [x,x*y,p]), x^2);
P3F2pxy!Evaluate(Gsing1g, [x,y,0]);
// "Old curve", now of genus 1:
Aff2F2 := AffineSpace(PF2uw);
Comp1 := Curve(Aff2F2, Evaluate(Gsing1g, [uu,ww,0]));
Genus(Comp1);
// Find its singular points:
SC1 := SingularSubscheme(Comp1);
Degree(SC1);
// Only one singular point.
Points(SC1);
// Tracing back, this is the other singularity on the original curve,
// which we will deal with later.
// The other component is:
P3F2pxy!Evaluate(Gsing1g, [0,y,p]);
// This is a smooth genus 1 curve (at least in this chart):
Comp2 := Curve(Aff2F2, Evaluate(Gsing1g, [0,ww,uu]));
Genus(Comp2);
IsNonsingular(Comp2);
//
// Let us determine the Euler factors for the two genus 1 curves.
PZ<T> := PolynomialRing(Integers());
// Note: ZetaFunction gives the zeta function
//       of the smooth projective model
Numerator(ZetaFunction(Comp1));
Numerator(ZetaFunction(Comp2));
// See Lemma 1.
//
// We need to check the intersection points of the two components
// if they are regular on the arithmetic surface.
P3F2pxy!Evaluate(Gsing1g, [0,y,0]);
// So the points are at x = p = 0, y = a  with  a^2 + a + 1 = 0.
// We look at the coefficients of x^0 p^0 with i+j <= 1:
Coefficient(Coefficient(Gsing1g, x, 0), p, 0);
// This is not in the ideal (x, p, y^2+y+1)^2, so the points are regular.
//
// We still have to check the other sub-chart
// for possible non-regular points.
// Second sub-chart: (x, y, p) <-- (y*x, y, p)
Gsing1h := ExactQuotient(Evaluate(Gsing1d, [y*x,y,p]), y^2);
P3F2pxy!Evaluate(Gsing1h, [x,y,0]);
// The "new" bit that was not visible in sub-chart 1 is where x = 0.
Evaluate($1, [0,yy,0]);
// So there is nothing in this chart of this component.
P3F2pxy!Evaluate(Gsing1h, [0,y,p]);
// Again, nothing.
// We also need to look at y = 0, since now 2 = p*x*y.
P3F2pxy!Evaluate(Gsing1h, [x,0,p]);
// This is the second genus 1 component again.
Evaluate($1, [0,yy,0]);
// Again, it does not meet this chart.
//
// Conclusion: The first singularity resolves into a new smooth curve
// of genus 1 intersecting the "old" curve (which is also of geometric
// genus 1) in two distinct points that are defined over F_4 and
// conjugate over F_2.
//
// Second singularity: <1,0> mod 2
// Shift to origin (mod 2)
Gsing2a := Evaluate(G, [U+1, W]);
[HomogeneousComponent(Gsing2a, d) : d in [0,1]];
// This shows that the point (u = 0, w = 1, 2 = 0) is non-regular,
// so we need to blow it up.
HomogeneousComponent(PF2uw!Gsing2a, 2);
// This has two distinct linear factors, so we have a node.
// The factors are defined over F_4 and conjugate.
//
// This already implies that the singularity will resolve into a chain
// of rational curves connecting two conjugate points on the resulting
// copy of the "old" curve. This is all we need to know here, since
// it implies that the dual graph of the special fiber has two independent
// loops which are both reversed by the action of Frobenius. 
//
// Together with the fact that there are two genus 1 curves in the special
// fiber, this shows that the connected component of the special fiber of
// the Neron model at p = 2 of the Jacobian is a semi-abelian variety,
// whose abelian part is a product of two elliptic curves, and whose toric
// part is a product of two copies of the norm-1 subgroup of F_4^*
// (at least up to isogeny). This gives the Euler factor
//   (1 + T)^2 (1 + T + 2T^2)^2   as given in Lemma 1; it also implies
// that the conductor exponent at p = 2 is 2.
//
// If one actually performs the necessary blow-ups, one sees that the
// chain consists of three rational curves.
//
// We now show that the Jacobian has trivial endomorphism ring (Lemma 2).
C := Curve(AffineSpace(P2Zuw), G);
EF5 := Numerator(ZetaFunction(BaseChange(C, GF(5)))); EF5;
EF7 := Numerator(ZetaFunction(BaseChange(C, GF(7)))); EF7;
// Note: ZetaFunction gives the zeta function
//       of the smooth projective model
PQtz<t,z> := PolynomialRing(Rationals(), 2);
// check criterion of Lemma 3
IsIrreducible(EF5); IsIrreducible(EF7);
Factorization(ExactQuotient(Resultant(Evaluate(EF5, z),
                                      Evaluate(EF5, t*z), z), (t-1)^8));
Factorization(ExactQuotient(Resultant(Evaluate(EF7, z),
                                      Evaluate(EF7, t*z), z), (t-1)^8));
// no cyclotomic factors
NF1 := NumberField(PZ!Reverse(Coefficients(EF5)));
NF2 := NumberField(PZ!Reverse(Coefficients(EF7)));
SF1 := Subfields(NF1); [Degree(s[1]) : s in SF1];
SF2 := Subfields(NF2); [Degree(s[1]) : s in SF2];
IsSubfield(SF1[1,1], NF2);
IsSubfield(SF2[1,1], NF1);
// ==> fields are linearly disjoint over Q
//
///////////////
// Section 3 //
///////////////
//
// Construct the curve in P^1 x P^1
Pr1 := ProjectiveSpace(Rationals(), 1);
Pr11<u1,u2,w1,w2> := DirectProduct(Pr1, Pr1);
// Homogenize the equation
Gpr11 := &+[MonomialCoefficient(G, U^i*W^j)*u1^i*u2^(3-i)*w1^j*w2^(3-j)
             : i, j in [0..3]];
Xpr11 := Curve(Pr11, Gpr11);
//
// Some functions don't work with curves in product projective spaces,
// so use Segre embedding to work in P^3.
Pr3<x00,x10,x01,x11> := ProjectiveSpace(Rationals(), 3);
segre := map<Pr11 -> Pr3 | [u1*w1,u2*w1,u1*w2,u2*w2]>;
Xpr3 := segre(Xpr11); Xpr3;
segreX := Restriction(segre, Xpr11, Xpr3);
//
// Look for rational points
ptsXpr3 := PointSearch(Xpr3, 100);
#ptsXpr3;
ptsXpr11 := [Points(pt @@ segreX) : pt in ptsXpr3];
// note: pulling back gives the preimage scheme, not a point.
assert forall{s : s in ptsXpr11 | #s eq 1};
ptsXpr11 := [s[1] : s in ptsXpr11]; ptsXpr11;
// These are points no. 8,7,3,5,0,4,9,2,6,1, in this order.
perm := [5, 10, 8, 3, 6, 4, 9, 2, 1, 7]; // ptsXpr11[perm[i+1]] = P_i
//
// Check their (t,c)-coordinates
cfnum := (-u1^3 - 2*u1^2*u2 + 5*u1*u2^2 - 10*u2^3)*u1*w1
           + (-u1^4 + 3*u1^3*u2 + 8*u1^2*u2^2 - 10*u1*u2^3 + 12*u2^4)*w2;
cfden := 4*u1^2*u2*(u1*w1+u1*w2-3*u2*w2);
cfunc := map<Xpr11 -> Pr1 | [cfnum, cfden]>;
// Check that c as given in the paper is correct
F2uw<u,w> := FunctionField(Rationals(), 2);
res := Evaluate(h6, [2/u-1, 0, 
                 Evaluate(cfnum, [u,1,w,1])/Evaluate(cfden, [u,1,w,1])]);
assert IsDivisibleBy(P2Zuw!Numerator(res), G);
// Now compute the c values; t = 2/u - 1 is clear.
Points(BaseScheme(cfunc)); // here the given expression is not defined
[Position(ptsXpr11, pt) : pt in $1];
[Position(perm, i) - 1 : i in $1]; // points P_0 and P_7
// Compute c for the others:
[cfunc(ptsXpr11[perm[i+1]]) : i in [0..9] | i notin $1];
// We have t(P_0) = infty, which implies c(P_0) = infty.
// For P_7, we have t(P_7) = -1; look at the polynomial relating t and c.
Roots(Evaluate(h6, [-1, 0, Polynomial(Rationals(), [0,1])]));
// So we must have c(P_7) = -2.
//
// The next task is to prove Lemma 4.
// We first compute the class groups of the reduction mod 3,5,7,11,13.
primes := [3, 5, 7, 11, 13]; 
Xred := [BaseChange(Xpr3, Bang(Rationals(), GF(p))) : p in primes];
Clgps := [<G, m> where G, m := ClassGroup(Xr) : Xr in Xred];
// Check that there is no rational torsion
Clnums := [&*[i : i in Invariants(pair[1]) | i ne 0] : pair in Clgps];
[Factorization(Clnums[3]), Factorization(Clnums[5])];
// This is for p = 7 and p = 13;
// 7^infty times the first and 13^infty times the second number
// are coprime.
//
// Set up the homomorphism Phi_S
// Get integral coordinates of P_0, ..., P_9
ptseqs := [Eltseq(ptsXpr3[perm[i+1]]) : i in [0..9]];
ptseqs[10] := [5*a : a in ptseqs[10]]; ptseqs;
ptseqs := [ChangeUniverse(s, Integers()) : s in ptseqs];
prodG, incls, projs := DirectProduct([pair[1] : pair in Clgps]);
Z10 := FreeAbelianGroup(10);
homs := [hom<Z10 -> Clgps[i,1] | [Divisor(Place(Xred[i]!pt)) @@ Clgps[i,2] 
                                   : pt in ptseqs]>  : i in [1..#Clgps]];
PhiS := hom<Z10 -> prodG | [&+[incls[i](homs[i](g)) : i in [1..#incls]]
                             : g in OrderedGenerators(Z10)]>;
Invariants(Image(PhiS));
// Since we know there is no torsion, this implies that the rank of
// our subgroup of J(Q) is at least 3 (or the rank of the Picard group
// of C at least 4).
Ker := Kernel(PhiS);
// Find small relations
L := Lattice(Matrix([Eltseq(Z10!k) : k in OrderedGenerators(Ker)]));
BL := Basis(LLL(L)); BL;
// We see six small elements in the kernel.
// Check that they correspond to principal divisors.
places := [Place(ptsXpr3[perm[i+1]]) : i in [0..9]];
divisors := [&+[BL[i,j]*places[j] : j in [1..#places]] : i in [1..6]];
assert forall{d : d in divisors | IsPrincipal(d)};
// The second return value of IsPrincipal is a function that has the
// given divisor, but we don't need these functions here.
//
// Set up the group G as a quotient of Z^10
Pic, mPic := quo<Z10 | [Z10!Eltseq(b) : b in BL[1..6]]>;
//
// We have now proved that the subgroup has at most rank 3,
// so the rank is 3.
// The relations in BL[1..6] are those given in the paper.
//
// Find the relation between the cusps:
sub<Z10 | [Eltseq(BL[i]) : i in [1..6]]>
  meet sub<Z10 | [Z10.i : i in [1..5]]>;
//
// Now approach Lemma 5.
// Take a basis of the subgroup of G
// that is in the kernel of reduction at 5.
hom5 := homs[2];
d1 := Z10![0, 0,0,0,0,0,0, 1, 0,-1]; // P_7 - P_9
d2 := Z10![1,-6,0,0,0,2,0, 1, 1, 1];
    // P_0 - 6 P_1 + 2 P_5 + P_7 + P_8 + P_9
d3 := Z10![1,-3,2,0,1,0,1,-1,-1, 0];
    // P_0 - 3 P_1 + 2 P_2 + P_4 + P_6 - P_7 - P_8
// Check that these divisors generate the kernel of reduction
assert mPic(sub<Z10 | d1,d2,d3>) eq mPic(Kernel(hom5));
//
// Turn them into divisors on the curve
Ds := [&+[s[i]*places[i] : i in [1..#places]] where s := Eltseq(d)
       : d in [d1,d2,d3]];
// The base point (as divisor)
base := Divisor(places[2]);
// Reduce D1, D2, D3 with respect to the base point
Dsred := [Reduction(D, base) : D in Ds];
[Degree(D) : D in Dsred];
// Dsred has the divisors D'_1, D'_2, D'_3.
// Find the characteristic polynomial of the u-values of the points 
// in these divisors:
Decomps := [Decomposition(D) : D in Dsred]; // [<place, multiplicity>, ...]
assert forall{D : D in Decomps | #D eq 1 and D[1,2] eq 1};
// ==> all D'_i are prime divisors of degree 4.
Dpts := [* RepresentativePoint(D[1,1]) : D in Decomps *];
// Extract u = x00/x10 and take minimal polynomial
charpols := [MinimalPolynomial(pt[1]/pt[2]) : pt in Dpts];
// Check that the geometric points in D'_i reduce to P'
// mod the prime above 5
[[Valuation(c, 5) : c in Coefficients(pol)] : pol in charpols];
// ==> except the last (which is the leading coefficient),
//     all coefficients are divisible by 5 ==> OK.
//
// Now find the differentials as power series in u, times du.
// We use z as the variable.
Pws<z> := LaurentSeriesAlgebra(Rationals());
// Express w = -1 + ... in terms of z
Roots(Evaluate(G, [z, Polynomial(Pws, [0,1])]));
// The first one of these has w = -1 mod z, so it is the correct one.
wseries := $1[1,1];
// Now the differentials are omega_0 = du/G_w etc.
omega0 := 1/Evaluate(Derivative(G, 2), [z, wseries]);
omega1 := z*omega0;
omega2 := wseries*omega0;
omega3 := z*omega2;
omegas := [omega0, omega1, omega2, omega3];
// Since P' is smooth on the reduction mod 5, these power series have
// 5-adically integral coefficients.
// We now integrate to obtain the logarithm series.
logs := [Integral(om) : om in omegas];
// When we plug in the power sums, the valuation of the term involving z^k
// will be at least ceil(k/4) - v_5(k).
// The terms up to k = 20 will therefore give a result up to O(5^5).
//
// Now compute the power sums. Note that
//  SUM_{k > 0} SUM_i u_i^k x^k = SUM_i u_i x/(1 - u_i x) = -x fr'(x)/fr(x)
// where  fr = PROD_i (1 - u_i x)  is the reciprocal of the polynomial
// that has the u_i as its roots.
powersums := [[Coefficient(s, i) : i in [1..20]]
                where s := -z*Derivative(f)/f
                where f := &+[Coefficient(pol, 4-j)*z^j : j in [0..4]]
                : pol in charpols];
// Plug into the logarithms:
Q5 := pAdicField(5);
mat := Matrix([[Q5!&+[Coefficient(l, i)*s[i] : i in [1..10]] + O(Q5!5^5)
                 : s in powersums] : l in logs])/5;
matker := KernelMatrix(ChangeRing(mat, GF(5))); matker;
// This shows that the reduction of the differential killing G is omega2.
//
// For dealing with (infty, 0), we need a little bit more precision.
matker2 := KernelMatrix(ChangeRing(mat, pAdicRing(5))); matker2;
// (over Q5, it doesn't work so well...)
//
// In the residue class of P_3, we take w as the uniformizer.
// Express the differentials in terms of w now. First u, which has a pole.
Roots(Evaluate(G, [Polynomial(Pws, [0,1]), z]));
// Only one root.
useries := $1[1,1];
omega0a := Derivative(useries)/Evaluate(Derivative(G, 2), [useries, z]);
omega1a := useries*omega0a;
omega2a := z*omega0a;
omega3a := z*omega1a;
omegasa := [omega0a, omega1a, omega2a, omega3a];
cofs := Coefficients(Integral(&+[(Rationals()!matker2[1,i])*omegasa[i]
                                  : i in [1..4]]));
cofs1 := [pAdicRing(5)!(cofs[i]*5^i) + O(pAdicRing(5)!5^4) : i in [1..5]];
[c/cofs1[1] : c in cofs1];
// ==> the power series that vanishes on C(Q) in that residue class is,
//     up to scaling,  l(tau) = tau (1 - 2*5 tau + O(5^2)) .
// In particular, its only root in Z_5 is tau = 0. (5*tau = w.)
//
///////////////
// Section 4 //
///////////////
//
// First, check the claims on the odd theta characteristics.
//
// Set up a scheme describing them
F<y0,y1,y2,z0,z1,z2,a,c,d> := PolynomialRing(Integers(), 9);
P2F<u,w> := PolynomialRing(F, 2);
g := Evaluate(G, [u,w]);
h1 := a + u + c*w + d*u*w; // a (1,1)-form
ru1 := Resultant(g, h1, w);
rw1 := Resultant(g, h1, u);
// both should be squares, up to scaling
du1 := ru1 - MonomialCoefficient(ru1,u^6)*(u^3+z2*u^2+z1*u+z0)^2;
dw1 := rw1 - MonomialCoefficient(rw1,w^6)*(w^3+y2*w^2+y1*w+y0)^2;
// Set up the ideal
eqns := Coefficients(du1) cat Coefficients(dw1);
I := ideal<F | eqns>;
// Verify that 1 + u + w is an odd thera characteristic
Variety(ChangeRing(ideal<F | I, a-1, c-1, d>, Rationals()));
//
// Consider scheme over F_5 and over F_13
I5 := ChangeRing(I, GF(5));
GB5 := GroebnerBasis(I5);
TotalDegree(GB5[#GB5]);
fact5 := Factorization(GB5[#GB5]);
[<TotalDegree(a[1]),a[2]> : a in fact5];
// ==> orbits under Fr_5 have lengths 1, 9 times 7, and 4 times 14
//
I13 := ChangeRing(I, GF(13));
GB13 := GroebnerBasis(I13);
TotalDegree(GB13[#GB13]);
fact13 := Factorization(GB13[#GB13]);
[<TotalDegree(a[1]),a[2]> : a in fact13];
// ==> orbits under Fr_13 have lengths 1 and 7 times 17
//
// If you are willing to wait for several minutes,
// you can also do it over Q:
// > time GB := GroebnerBasis(ChangeRing(I, Rationals()));
// > TotalDegree(GB[#GB]);
// > fact := Factorization(GB[#GB]);
// > [<TotalDegree(a[1]),a[2]> : a in fact];
//
// Set up the symplectic group
Sp := SymplecticGroup(8, 2);
// We know that Fr_5 and Fr_13 are elements of order 14 and 17, resp.
maxSp := [M : M in MaximalSubgroups(Sp) | IsDivisibleBy(M`order, 14*17)];
assert #maxSp eq 1; // only one such subgroup
// Extract the subgroup
Gamma := maxSp[1]`subgroup;
#Sp/#Gamma;
// Since the full symplectic group acts transitively on J[2] \ {0}, but
// we have an orbit of size 119, the Galois group must map into Gamma.
// Check that it cannot be smaller.
maxGamma := MaximalSubgroups(Gamma);
maxG1 := [M : M in maxGamma | IsDivisibleBy(M`order, 14*17)];
assert #maxG1 eq 1;
Gamma1 := maxG1[1]`subgroup;
maxGamma1 := MaximalSubgroups(Gamma1);
assert IsEmpty([M : M in maxGamma1 | IsDivisibleBy(M`order, 14*17)]);
// Check that Gamma1 is not OK,
// since it does not have elements of order 14.
assert forall{cl : cl in ConjugacyClasses(Gamma1) | Order(cl[3]) ne 14};
// So the image of Galois must be Gamma.
// Find the orbits on J[2].
orbs := Orbits(Gamma);
[#o : o in orbs];
// ==> orbits of lengths 1, 119, 136
//
// Find subgroups of smallest index in Gamma:
[#Gamma/M`order : M in maxGamma];
// ==> smallest are 2 and 119; check subgroups of Gamma1
//     (with (Gamma:Gamma1)=2)
[#Gamma/M`order : M in maxGamma1];
// ==> nothing smaller than before; the index-2 subgroup does not give
//     a faithful permutation representation
// ==> need at least degree 119.
//
// Now proceed to the L-series computations.
cond := 2^2*bigp; // the conductor
// The precision we will use. Increase to get higher precision.
// (The computation will then take longer.)
prec := 5; 
Ls := LSeries(2, [0,0,0,0,1,1,1,1], cond, 0
               : Sign := -1, Precision := prec);
numcofs := LCfRequired(Ls); numcofs;
// For a later check that the sign is correct:
Ls1 := LSeries(2, [0,0,0,0,1,1,1,1], cond, 0
                : Sign := 1, Precision := prec);
//
// Now set up a list of Euler factors up to the required precision.
// We know the factor at 2; it is (1 + T)^2 (1 + T + 2 T^2)^2 .
// All other relevant primes are good.
pr4 := [p : p in [3..Floor(numcofs^(1/4)) by 2] | IsPrime(p)];
pr3 := [p : p in [pr4[#pr4]+2..Floor(numcofs^(1/3)) by 2] | IsPrime(p)];
pr2 := [p : p in [pr3[#pr3]+2..Floor(numcofs^(1/2)) by 2] | IsPrime(p)];
pr1 := [p : p in [pr2[#pr2]+2..numcofs by 2] | IsPrime(p)];
// For the smallest primes, use the ZetaFunction to get the Euler factor.
XZ := BaseChange(Xpr3, Bang(Rationals(), Integers()));
EFs4 := [<2, (1+T)^2*(1+T+2*T^2)^2>]
          cat [<p, Numerator(ZetaFunction(BaseChange(XZ, GF(p))))>
                : p in pr4];
EFs3 := [<p, Numerator(ZetaFunction(BaseChange(XZ, GF(p))))> : p in pr3];
// For reasons of efficiency, we use an affine plane curve.
// Note that the five points at infinity are all rational and remain 
// distinct mod p for all odd p.
XZaff := Curve(AffineSpace(Parent(G)), G); 
// For the larger primes, count points in C(F_p), plus perhaps in C(F_p^2).
ptcounts2 := [<p, [#Points(BaseChange(XZaff, GF(p)), GF(p,n))+5
               : n in [1..2]]> : p in pr2];
EFs2 := [<a[1], 1 + (a[2,1]-a[1]-1)*T
            + (ExactQuotient(a[2,2]+a[2,1]^2,2)-a[2,1]*(a[1]+1)+a[1])*T^2>
           : a in ptcounts2];
EFs1 := [<p, 1+(#Points(BaseChange(XZaff,GF(p)),GF(p))-p+4)*T> : p in pr1];
// Set up a map
EFmap := pmap<Integers() -> PZ | EFs4 cat EFs3 cat EFs2 cat EFs1>;
// Set coefficient function
LSetCoefficients(Ls, func<p, d | EFmap(p)>);
LSetCoefficients(Ls1, func<p, d | EFmap(p)>);
// Now check the functional equation
CheckFunctionalEquation(Ls);
// This is reasonably close to zero relative to the given precision.
CheckFunctionalEquation(Ls1);
// This is clearly nonzero.
// ==> the negative sign is correct.
//
// Now compute L(1), L'(1), L''(1), L'''(1)
Evaluate(Ls, 1);
Evaluate(Ls, 1 : Derivative := 1, Leading);
Evaluate(Ls, 1 : Derivative := 2, Leading);
Evaluate(Ls, 1 : Derivative := 3, Leading);
// The first three are zero (up to the given precision),
// the last clearly isn't
// ==> ord_{s=1} L(C, s) = 3.