(* This file is organized according to the section 
structure of the paper, and is intended to be human-readable.
The other files define:

p40:  a degree 40 polynomial in Q[a,b,c,d,x] with 1673 terms

mat240 and matstar240: two 240-by-240 matrices with entries in Q[a,b,c,d,s,t,u,v].
Together they require about 5 megabytes to store.  

ABCD: a list of four polynomials in Q[a,b,c,d,s,t,u,v] which together take
145 megabytes to store

ABCDstar: four similar polynomials where (s,t,u,v) = (1,0,0,0) and 
so only require 18 kilobytes to store.  

Thus the file ABCD, which corresponds to the main theorem of the 
paper, is the bulkiest among the ancillary files. *)

(* As a sample use, for the genus two curve 

X: y^2 = x^5+x^3+2 x^2 + 3 x + 4,

the quintic defining the full parametric family of genus two curves

Y: y^2 = x^5+ A x^3+B x^2 + C x + D,

with a symplectic isomorphism Jac(X)[3] -> Jac(Y)[3] is obtained by
typing 

quintic[ABCD[{1,2,3,4,s,t,u,v}]].

One member of this family is given by 

quintic[ABCD[{1,2,3,4,5,6,7,8}]].

Because there are so many terms in the generic ABCD[{a,b,c,d,s,t,u,v}], 
each of these last commands takes 
a little over ten seconds to run.  A quicker version in the
second case would be

quintic[ABCDmat[{1,2,3,4,5,6,7,8}]].

The commands ABCD and ABCDmat come from Section 3.5, and more
about them is given in the comments there.
*)

(*################################################################*)
 
(* Section ONEone.  Overview. *)

quintic[{a_,b_,c_,d_}] := x^5 + a*x^3 + b*x^2 + c*x + d


(*################################################################*)

(* Section TWOone.  Elliptic curves and their 3-torsion. *)

Delta[{a_,b_}] := 2^4(-4*a^3 - 27*b^2)

F[{a_,b_},z_] := z^8 + 18*a*z^4 + 108*b*z^2 - 27*a^2


(*################################################################*)

(* Section TWOthree.  Rings of Invariants.   *)

rels23= {w,a,b}-{u^3/3+u^2*z + u*z^2,w*z/9,(w^2-6*w*z^3 - 3*z^6)/324}

check23 := Factor[Resultant[rels23[[2]],rels23[[3]],w]/F[{a,b},z]] 
(* returns a constant, confirming proportionality *)


(*################################################################*)

(* Section TWOfour.  Covariants and contravariants.   *)

alphas[{a_,b_}] := {z, (w + z^3)/6} (* returns alpha1,alpha3 *)

betas[{a_,b_}] := {(w - z^3)/2, (5*w*z^2 + 3*z^5)/18} (* returns beta3,beta5 *)


(*################################################################*)

(* Section TWOfive.  New coefficients.   *)

AB[{a_,b_,s_,t_}] := 
{(3*a*s^4 + 18*b*s^3*t - 6*a^2*s^2*t^2 - 6*a*b*s*t^3 - a^3*t^4 - 9*b^2*t^4)/3,
(9*b*s^6 - 12*a^2*s^5*t - 45*a*b*s^4*t^2 - 90*b^2*s^3*t^3 + 15*a^2*b*s^2*t^4 - 
  4*a^4*s*t^5 - 18*a*b^2*s*t^5 - 3*a^3*b*t^6 - 18*b^3*t^6)/9}

ABstar[{a_,b_,s_,t_}] := 
{108*a^3*s^4 - 243*b^2*s^4 - 864*a^2*b*s^3*t - 72*a^4*s^2*t^2 + 
  1458*a*b^2*s^2*t^2 + 144*a^3*b*s*t^3 - 972*b^3*s*t^3 - 4*a^5*t^4 - 
  135*a^2*b^2*t^4
,
-2*(972*a^3*b*s^6 + 729*b^3*s^6 + 432*a^5*s^5*t - 
   4860*a^2*b^2*s^5*t - 2700*a^4*b*s^4*t^2 + 10935*a*b^3*s^4*t^2 + 
   7560*a^3*b^2*s^3*t^3 - 7290*b^4*s^3*t^3 + 180*a^5*b*s^2*t^4 - 
   8505*a^2*b^3*s^2*t^4 + 16*a^7*s*t^5 - 108*a^4*b^2*s*t^5 + 4374*a*b^4*s*t^5 - 
   20*a^6*b*t^6 - 135*a^3*b^3*t^6 - 1458*b^5*t^6)
}

check25 := 
Factor[Resultant[Numerator[Factor[(u - {s,t}.alphas[{a,b}]) /. (w->9*a/z)]],
F[{a,b},z],z]/
F[AB[{a,b,s,t}],u]] (* returns a quantity independent of u *)


check25star := 
Factor[Resultant[Numerator[Factor[(u - {s,t}.betas[{a,b}] ) /. (w->9*a/z)]],
F[{a,b},z],z]/
F[ABstar[{a,b,s,t}],u]] (* returns a quantity independent of u *)


(*################################################################*)

(* Section TWOsix.   Geometric summary.   *)

check26 := Factor[Delta[AB[{a,b,s,t}]]/Delta[{a,b}]] 
(* a perfect cube *)

check26star := 
Factor[Delta[ABstar[{a,b,s,t}]]*Delta[{a,b}]] 
(* a perfect cube *)


(*################################################################*)

(* Section TWOseven.  Finding (s,t).  *)

check27a := {-27,-162}-ABstar[{-1,0,-1/2,3/2}] 
(* returns {0,0}, confirming the proper identification 
of an antisymplectic isomorphism
between the 3-torsion of two elliptic curves *)

check27b := {5805,-285714} - ABstar[{5805,-285714,435/254016,11/254016}]
(* returns {0,0}, confirming the proper identification of an 
antisymplectic automorphism
on the 3-torsion of X_0(14) *)


(*################################################################*)

(* Section THREEone.  Weierstrass curves and their 3-torsion. *)

Delta[{a_,b_,c_,d_}] := 2^8*Discriminant[quintic[{a,b,c,d}],x]

(* p40 defined in the file p40.txt *)

F[{a_,b_,c_,d_},z_] := p40[{a,b,c,d},z^6]

check31 := Factor[F[{a,b,c,d},0]/(3^12*Delta[{a,b,c,d}]^2)] (* returns1 *)


(*################################################################*)

(* Section THREEthree.  Rings of Invariants.   *)

rels33a = {p,q,r}-
{z1^6 + z2^6 + z3^6 - 10*(z1^3*z2^3 + z1^3*z3^3 + z2^3*z3^3), (z1^3 - z2^3)*(z2^3 - z3^3)*(z3^3 - z1^3), 
 (z1^3 + z2^3 + z3^3)*(216*z1^3*z2^3*z3^3 + (z1^3 + z2^3 + z3^3)^3)}

rels33b = ({2^4*3^7*5*a , 2^6*3^9*5^2*b  , 2^8*3^12*5^3*c  , 2^10*3^16*5^5*d}
-
{-p^2 - 5*r + 1320*q*z^3 - 132*p*z^6 - 6*z^12, p^3 - 400*q^2 - 5*p*r - 680*p*q*z^3 + 
  (323*p^2 - 255*r)*z^6 - 7480*q*z^9 + 68*p*z^12 - 4*z^18, 2*p^4 - 800*p*q^2 - 5*p^2*r + 
  (320*p^2*q - 3000*q*r)*z^3 + (-722*p^3 + 175200*q^2 + 990*p*r)*z^6 + 33040*p*q*z^9 + 
  (-953*p^2 + 3495*r)*z^12 + 15720*q*z^15 + 268*p*z^18 - 3*z^24, 
 13*p^5 - 6000*p^2*q^2 - 25*p^3*r + (21600*p^3*q - 9600000*q^3 - 45000*p*q*r)*z^3 + 
  (11790*p^4 - 4572000*p*q^2 - 37575*p^2*r + 28125*r^2)*z^6 - (247200*p^2*q + 945000*q*r)*z^9 + 
  (37155*p^3 + 234000*q^2 - 150075*p*r)*z^12 - 214200*p*q*z^15 + (30855*p^2 - 143775*r)*z^18 + 
  354600*q*z^21 + 2340*p*z^24 - 12*z^30})


(*################################################################*)

(* Section THREEfour.  Covariants and contravariants.   *)

alphas[{z_,p_,q_,r_}] := (* returns (alpha1,alpha7,alpha13,alpha19) *)
{z, -(z*(-7*p + 3*z^6))/540, (z*(-3*p^2 + 11*r + 216*q*z^3 + 72*p*z^6))/11664,
 -(z*(-p^3 + 468*q^2 + p*r + 24*p*q*z^3 + 6*p^2*z^6 - 66*r*z^6 + 288*q*z^9 + 12*p*z^12))/944784}

betas[{z_,p_,q_,r_}] := (* returns (beta11,beta17,beta23,beta29) *)
{
-(z^2*(-55*q + 11*p*z^3 + z^9))/2430,
 -(z^2*(340*p*q - 323*p^2*z^3 + 255*r*z^3 + 11220*q*z^6 - 136*p*z^9 + 12*z^15))/5248800,
 -(z^2*(-80*p^2*q + 750*q*r + 361*p^3*z^3 - 87600*q^2*z^3 - 495*p*r*z^3 - 
24780*p*q*z^6 + 953*p^2*z^9 - 3495*r*z^9 - 19650*q*z^12 -
     402*p*z^15 + 6*z^21))/1417176000,
 -(z^2*(-720*p^3*q + 320000*q^3 + 1500*p*q*r - 786*p^4*z^3 + 304800*p*q^2*z^3 +
 2505*p^2*r*z^3 - 1875*r^2*z^3 + 24720*p^2*q*z^6 +
     94500*q*r*z^6 - 4954*p^3*z^9 - 31200*q^2*z^9 + 20010*p*r*z^9 + 35700*p*q*z^12 - 
6171*p^2*z^15 + 28755*r*z^15 - 82740*q*z^18 -
     624*p*z^21 + 4*z^27))/1530550080000}


(*################################################################*)

(* Section THREEfive.  New coefficients. *)

(* We present three programs which in principle agree on inputs
{a,b,c,d,s,t,u,v} in Q^8:

ABCDnum[{a,b,c,d,s,t,u,v}] = ABCDmat[{a,b,c,d,s,t,u,v}] = ABCD[{a,b,c,d,s,t,u,v}]

There are also starred versions 

ABCDstarnum[{a,b,c,d,s,t,u,v}] = ABCDstarmat[{a,b,c,d,s,t,u,v}] = ABCDstar[{a,b,c,d,s,t,u,v}],

although ABCDstar works only if {s,t,u,v}={1,0,0,0} *)

(* The numerical approach is not described explicitly in the paper, although its main ingredients are. 
It has the advantage that the code is short, as it does not need any of the data
outside of this file.   
Indeed it uses only the alphas and betas from S3.4, and the commands below through ABCDstarnum.  
It requires a, b, c, and d to be rational numbers, 
but allows s, t, u, v to be free parameters.

As an example of timing, ABCDnum[{1,2,3,4,5,6,7,8}] takes 35 seconds to run. 
After that, the numerics for any
ABCDnum[{1,2,3,4,s,t,u,v}] have been stored, decreasing run times.  For inputs

{1,2,3,4,5,6,7,9},
{1,2,3,4,5,6,7,v},
{1,2,3,4,5,6,u,v},
{1,2,3,4,5,t,u,v},
{1,2,3,4,s,t,u,v},

the run times for ABCDnum are about 0.2, 1, 13, 60, and 62 seconds.  Run times for ABCDstarnum are
very similar. 
*)

prec=500 (* this is more than enough for the examples 
in the paper and file, but would of course need to be increased if
inputs have sufficiently large heights *)

rootvecs80[{a_,b_,c_,d_}] := rootvecs80[{a,b,c,d}] =
{y,p,q,r} /.
NSolve[
{174960*a + p^2 + 5*r - 1320*q*y + 132*p*y^2 + 6*y^4 == 0, 31492800*b - p^3 + 400*q^2 + 5*p*r + 680*p*q*y - (323*p^2 - 255*r)*y^2 +
   7480*q*y^3 - 68*p*y^4 + 4*y^6 == 0, 17006112000*c - 2*p^4 + 800*p*q^2 + 5*p^2*r - (320*p^2*q - 3000*q*r)*y -
   (-722*p^3 + 175200*q^2 + 990*p*r)*y^2 - 33040*p*q*y^3 - (-953*p^2 + 3495*r)*y^4 - 15720*q*y^5 - 268*p*y^6 + 3*y^8 == 0,
 137749507200000*d - 13*p^5 + 6000*p^2*q^2 + 25*p^3*r - (21600*p^3*q - 9600000*q^3 - 45000*p*q*r)*y -
   (11790*p^4 - 4572000*p*q^2 - 37575*p^2*r + 28125*r^2)*y^2 + (247200*p^2*q + 945000*q*r)*y^3 -
   (37155*p^3 + 234000*q^2 - 150075*p*r)*y^4 + 214200*p*q*y^5 - (30855*p^2 - 143775*r)*y^6 - 354600*q*y^7 - 2340*p*y^8 + 12*y^10 == 0}
,{y,p,q,r},prec]

rootvecs240[{a_,b_,c_,d_}] := Flatten[Table[{r[[1]]^(1/3),r[[2]],r[[3]],r[[4]]}*{Exp[2 Pi I j/3],1,1,1},
{r,rootvecs80[{a,b,c,d}]},{j,0,2}],1]

ABCDfromtraces[{t2_,t3_,t4_,t5_}] :=
Rationalize[Expand[{-t2/30240, -t3/7862400, (3667*t2^2 - 5600*t4)/9390915072000, (2521*t2*t3 - 2688*t5)/886312627200000}]]

ABCDnum[{a_,b_,c_,d_,s_,t_,u_,v_}] :=
Module[{covariants,roots},
covariants = Map[alphas,rootvecs240[{a,b,c,d}]];
roots = Table[cov.{s,t,u,v},{cov,covariants}];
ABCDfromtraces[{Total[roots^12],Total[roots^18],Total[roots^24],Total[roots^30]}/6]
]

ABCDstarnum[{a_,b_,c_,d_,s_,t_,u_,v_}] :=
Module[{contravariants,roots},
contravariants = Map[betas,rootvecs240[{a,b,c,d}]];
roots = Table[contra.{s,t,u,v},{contra,contravariants}];
ABCDfromtraces[{Total[roots^12],Total[roots^18],Total[roots^24],Total[roots^30]}/6]
]


(* The approach using 240-by-240 matrices is presented in the paper.  The paper describes how to 
construct the basic matrices M and M^*.  This file does not implement the construction. Instead, 
the matrices can be found in the file mat240-and-matstar240.txt.
  
The matrix approach gives proved answers, there being no precision issues.  However run times 
are much slower than the those of the numerical approach.  
For example ABCDmat[{1,2,3,4,5,6,7,8}] takes 2.4 seconds while ABCDmat[{1,2,3,4,5,6,7,v}] takes
152 seconds.  A key advantage is now that the first four arguments are allowed to be 
unspecified.  For example, ABCDmat[{1,2,3,d,5,6,7,8}] takes 87 seconds.  For ABCDstarmat the
timings for the same inputs are 1.8, 183, and 150 seconds.  
These very slow run times are why the decomposition
of M^k according to the multinomial theorem is highlighted in the text.
*) 

(* mat240 and matstar240 are defined the file mat240-and-matstar240.txt. *)

trace[mat1_,mat2_] := Module[{tmat2},tmat2 = Transpose[mat2];
Sum[Expand[mat1[[i]].tmat2[[i]]],{i,1,Length[mat1]}]]

ABCDmat[{a_,b_,c_,d_,s_,t_,u_,v_}] := Module[{mat,mat3,mat6,mat9,mat12,mat15},
mat = mat240[{a,b,c,d,s,t,u,v}];
mat3 = Expand[Expand[mat.mat].mat];
mat6 = Expand[mat3.mat3];
mat9 = Expand[mat3.mat6];
mat12 = Expand[mat6.mat6];
mat15 = Expand[mat12.mat3];
ABCDfromtraces[
{
trace[mat6,mat6],
trace[mat9,mat9],
trace[mat12,mat12],
trace[mat15,mat15]
}/6]]

ABCDstarmat[{a_,b_,c_,d_,s_,t_,u_,v_}] := Module[{mat,mat3,mat6,mat9,mat12,mat15},
mat = Expand[matstar240[{a,b,c,d,s,t,u,v}]];
mat3 = Expand[Expand[mat.mat].mat];
mat6 = Expand[mat3.mat3];
mat9 = Expand[mat3.mat6];
mat12 = Expand[mat6.mat6];
mat15 = Expand[mat12.mat3];
ABCDfromtraces[
{
trace[mat6,mat6],
trace[mat9,mat9],
trace[mat12,mat12],
trace[mat15,mat15]
}/6]]

(* The command ABCD just gives the stored answer. 
 There are retrieval issues with so many terms,
and so completely specialized cases like ABCD[{1,2,3,4,5,6,7,8}] 
take a long time, namely about 14 seconds.  
But the fully generic case ABCD[{a,b,c,d,s,t,u,v}] also takes about 14 seconds.  *)

(* ABCD and ABCDstar are in the files ABCD.txt and ABCDstar.txt respectively. *)


(*################################################################*)

(* Section THREEseven.  Finding (s,t,u,v). *)

(* The material in this file correspoding to 3.7 is long but completely self-contained, 
except that uses quintic and Delta and requires a precision prec to be set; the default given earlier is 500.  
The main command is findisos, which 
appears halfway through the discussion after the introduction
of some other commands.   findisos inputs coefficient vectors
{a,b,c,d} and {aa,bb,cc,dd} in Q^4 defining Weierstrass curves X and Y 
respectively.  As a preliminary, it gets two new coefficient
vectors {aa1,bb1,cc1,dd1} and {aa2,bb2,cc2,dd2} defining Y1
and Y2 isomorphic to Y. Here {aai,bbi,cci,ddi} = {aa,bb,cc,dd}*ri^{4,6,8,10} 
where the ri are rescaling factors.  These rescaling factors 
are chosen so that for the  adjusted curves
Delta_X/Delta_Y1 and Delta_X Delta_Y2 are rational cubes. 

findisos returns {{r1,stuv1list},{r2,stuv2list}}. 
Each element of stuvilist is a nonzero element {s,t,u,v} of Q^4. 
The stuvilists contain up to eight vectors.  Each list
is stable under multiplication by -1 and so is even
in length.  If {a,b,c,d} and {aa,bb,cc,dd} have
nonisomorphic 3-torsion then both lists have length
0.  Otherwise, at least one list has positive length,  
with i=1 corresponding to symplectic isomorphisms
and i=2 corresponding to antisymplectic isomorphisms.
The generic behavior in the setting of isomorphic
3-torsion is that one list has length 0 and the other
length 2.  

The discussion concludes with examples where there are indeed
isomorphisms on 3-torsion.  Typing e.g. 
examp1 or examp2[17] will
return output from findisos after about
ten seconds.  
*)

(* eqs40 is called only with arguments {a,b,c,d}.  
It basically sets up the equations expressing the the G-invariants 
{a,b,c,d} in terms of the H-invariants {z,p,q,r} 
Ultimately there are 240 solutions, but they come in packets of 6, so
we get the 40 packets first, using quantities (y,u6,v6,u12) closely
related to (z,p,q,r).
*)

eqs40[{a_, b_, c_, d_}] := 
{a + (20*u12 + u6^2 + 22*u6*y - 220*v6*y + y^2)/29160, 
 b + (30*u12*u6 + u6^3 + 1530*u12*y - 17*u6^2*y + 170*u6*v6*y + 100*v6^2*y - 
    17*u6*y^2 + 1870*v6*y^2 + y^3)/7873200, 
 c + (120*u12*u6^2 + 3*u6^4 - 23760*u12*u6*y - 268*u6^3*y + 72000*u12*v6*y + 
    2680*u6^2*v6*y + 800*u6*v6^2*y - 83880*u12*y^2 - 2542*u6^2*y^2 - 
    33040*u6*v6*y^2 - 175200*v6^2*y^2 - 268*u6*y^3 - 15720*v6*y^3 + 3*y^4)/
   17006112000, d + (50*u12*u6^3 + u6^5 - 1350000*u12^2*y - 37350*u12*u6^2*y - 
    195*u6^4*y + 90000*u12*u6*v6*y + 1950*u6^3*v6*y + 500*u6^2*v6^2*y + 
    300150*u12*u6*y^2 + 9410*u6^3*y^2 + 1890000*u12*v6*y^2 + 
    99350*u6^2*v6*y^2 + 381000*u6*v6^2*y^2 + 800000*v6^3*y^2 + 287550*u12*y^3 + 
    9410*u6^2*y^3 + 17850*u6*v6*y^3 - 19500*v6^2*y^3 - 195*u6*y^4 - 
    29550*v6*y^4 + y^5)/11479125600000}

(* variants240 is only called with the 
arguments {a,b,c,d}.  It returns a list of 240 vectors
in C^8.   On each list, eight of the
vectors are in R^8.  variantsR picks 
the first of these real vectors, which
we call {alpha1,alpha7,alpha13,alpha19,beta11,beta17,beta23,beta29}
*)

variants240[{a_, b_, c_, d_}] :=
  Module[{list40, list240},
   list40 = {y, u6, v6, u12} /.
     NSolve[Table[eq == 0, {eq, eqs40[{a, b, c, d}]}], {y, u6, v6,
       u12}, prec];
   list240 =
    Flatten[Table[{Exp[(2 Pi I j)/6] u[[1]]^(1/6), u[[2]],
       u[[3]] Exp[(2 Pi I j)/2] Sqrt[u[[1]]], u[[4]]}, {u,
       list40}, {j, 0, 5}], 1];
   Table[
{x, (x*(7*u6 - 3*x^6))/540, (x*(33*u12 + u6^2 + 27*u9*x^3 + 9*u6*x^6))/1458, 
 -(x*(2*u12*u6 + 39*u9^2 + 2*u6*u9*x^3 - 132*u12*x^6 - 5*u6^2*x^6 + 24*u9*x^9 + 
     u6*x^12))/78732, (x^2*(55*u9 - 11*u6*x^3 - x^9))/2430, 
 (x^2*(-85*u6*u9 - 1530*u12*x^3 + 17*u6^2*x^3 - 2805*u9*x^6 + 34*u6*x^9 - 
    3*x^15))/1312200, 
 -(x^2*(9000*u12*u9 + 335*u6^2*u9 - 5940*u12*u6*x^3 - 67*u6^3*x^3 - 
     43800*u9^2*x^3 - 12390*u6*u9*x^6 - 41940*u12*x^9 - 1271*u6^2*x^9 - 
     9825*u9*x^12 - 201*u6*x^15 + 3*x^21))/708588000, 
 (x^2*(-9000*u12*u6*u9 - 195*u6^3*u9 - 80000*u9^3 + 270000*u12^2*x^3 + 
    7470*u12*u6^2*x^3 + 39*u6^4*x^3 - 76200*u6*u9^2*x^3 - 567000*u12*u9*x^6 - 
    29805*u6^2*u9*x^6 - 120060*u12*u6*x^9 - 3764*u6^3*x^9 + 7800*u9^2*x^9 - 
    8925*u6*u9*x^12 - 172530*u12*x^15 - 5646*u6^2*x^15 + 20685*u9*x^18 + 
    156*u6*x^21 - x^27))/382637520000}
/. {x -> li[[1]], u6 -> li[[2]], u9 -> li[[3]],
      u12 -> li[[4]]}, {li, list240}]]

variantsR[{a_,b_,c_,d_}] := Select[variants240[{a,b,c,d}],Im[Chop[N[#[[1]],prec/2]]]==0 &][[1]]

(* newroots240 is only called with the coefficient vectors 
{aai,bbi,cci,ddi} and gives the first column of 
variants240[{aai,bbi,cci,ddi}].  newroots8R picks
out the eight which are real, called Z1,...,Z8 in 
the i=1 case and Zstar1,...,Zstar8 in the i=2 case *)

newroots240[{aa_,bb_,cc_,dd_}] := 
  Module[{(* list40 *)},
   list40 = y /.  NSolve[Table[eq == 0, {eq, eqs40[{aa, bb, cc, dd}]}], {y, u6, v6, u12}, prec];
    Flatten[Table[Exp[(2 Pi I j)/6] u^(1/6),{u,list40}, {j, 0, 5}], 1]]

newroots8R[{aa_,bb_,cc_,dd_}] := Select[newroots240[{aa,bb,cc,dd}],Chop[Im[N[#,prec/2]]]==0 &]

(* scale factors *)

cubefree[n_] := Module[{fi},fi = FactorInteger[Abs[n]];
Product[f[[1]]^Mod[f[[2]],3],{f,fi}]]

r1[{abcd_,aabbccdd_}] :=  cubefree[Delta[abcd]/Delta[aabbccdd]]

r2[{abcd_,aabbccdd_}] :=  cubefree[1/Delta[abcd]/Delta[aabbccdd]]

(* findrel uses Mathematica's LatticeReduce to find 
an integer vector w with modest height entries so that
w.v is near 0 for the inputed real vector.  findrel2 
reformats the answer slightly. *)

findrel[v_] := Module[{hold, m, tab}, hold = Round[10^prec v];
  m = IdentityMatrix[Length[hold]];
  tab = Table[Join[m[[i]], {hold[[i]]}], {i, 1, Length[hold]}];
  Drop[LatticeReduce[tab][[1]], -1]]

reformat[vec_] := -Drop[vec,-1]/vec[[-1]]

findrel2[v_] := reformat[findrel[v]]

(* findisos is explained at the top of this subsection.  In practice,
meaningless results have height about 10^(prec/4) and meaningful
ones have heights much less than 10^(prec/8), so that
the identification of correct solutions seems robust. *)

adjust[lis_,z_] := lis*(z^{4,6,8,10})

height[r_] := Max[Abs[Numerator[r]],Denominator[r]]

findisos[{a_,b_,c_,d_},{aa_,bb_,cc_,dd_}] := Module[{invs,nr1,nr2},
invs = variantsR[{a,b,c,d}];
nr1 = newroots8R[adjust[{aa,bb,cc,dd},r1[{{a,b,c,d},{aa,bb,cc,dd}}]]];
nr2 = newroots8R[adjust[{aa,bb,cc,dd},r2[{{a,b,c,d},{aa,bb,cc,dd}}]]];
cands1 = Table[Join[Take[invs,4],{nr}],{nr,nr1}];
cands2 = Table[Join[Take[invs,-4],{nr}],{nr,nr2}];
{{r1[{{a,b,c,d},{aa,bb,cc,dd}}],
Select[Map[findrel2,cands1],height[#[[1]]]<10^(prec/8) &]},
{r2[{{a,b,c,d},{aa,bb,cc,dd}}],Select[Map[findrel2,cands2],height[#[[2]]]<10^(prec/8) &]}}]

(* EXAMPLE 1: an atypical case.  By setting (z1,z2,z3,z) = (1,2,4,10) and 
tracing through the equations of Section 3 one gets
the coefficient list split124x.  The degree 240 polynomial
F[split124x,z] is separable with eight linear factors and 116 quadratic 
factors splitting over E=Q(sqrt(-3)).  In other words, 
the always present situation for 
C/R is here seen at the level of E/Q.    
examp1 gives eight symplectic isomorphisms 
(two of which are just {1,0,0,0} and {-1,0,0,0})
and eight antisymplectic isomorphisms.
This example is analogous to
the example of X_0(14) in the paper's S2.7.
*)

split124x = 
{-968525982721/29160, 
-984465932430343681/7873200, 
-1237329868861009891031041/5668704000, 
 -2486798498191636627695169804801/11479125600000}

examp1 := findisos[split124x,split124x]

(* EXAMPLE 2:  39 individual examples with small minimal discriminant coming 
from the LMFDB.   

As of February 2020 the LMFDB has a presumably near-complete list 
of all genus two curves with minimal absolute discriminant <=1000000. 
For exactly 39 Galois representation
rhobar with full image GSp_4(F_3) are there at least two Weierstrass curves
having non-isogenous Jacobians with 3-torsion given by rhobar.  
Always there are only two isogeny classes. 
Always the two isogeny classes have different conductor.  With four exceptions, 
all classes under consideration contain exactly one Weierstrass curve. 
These classes are the one with lower conductor in cases 1, 6, and 8 
and the one with larger conductor in case 30 on the lmfdb list. 
All exceptional classes contain exactly two Weierstrass curves.

Entry i on the smallcond gives two or exceptionally three pairs, sorted 
by conductor.  The pairs are in the form {{N,letter,Delta,number},{a,b,c,d}}.
Here N.letter is the LMFDB label for the isogeny class and N.letter.Delta.number
is the LMFDB label for the curve.  In the four cases where an isogeny class
is represented by two curves, 3-torsion representations are symplectically
isomorphic.  When one changes conductor the 3-torsion representations
are symplectically isomorphic in 34 cases.  In cases 2, 16, 18, 28, and 34
they are antisymplectically isomorphic.  

To see findisos applied to case i, ignoring the middle curve if there 
are three listed curves, type examp2[i] *)

examp2[i_] := findisos[smallcond[[i,1,2]],smallcond[[i,-1,2]]]

smallcond = 
{{{{461, a, 461, 1}, {14000, 7320000, 705600000, 20010400000}}, 
  {{461, a, 461, 2}, {-394000, -13560000, 27129600000, -2512632800000}}, 
  {{1844, b, 29504, 1}, {24750, -3502500, -5634375, 7406237500}}}, 
 {{{644, b, 14812, 1}, {54750, -552500, -53384375, 847487500}}, 
  {{2576, a, 288512, 1}, {-66000, 1720000, 1141600000, -61209600000}}}, 
 {{{976, a, 999424, 1}, {-152250, 54442500, -15637009375, 2299537412500}}, 
  {{1952, b, 999424, 1}, {44000, -22080000, 509600000, 198630400000}}}, 
 {{{1051, a, 1051, 1}, {-4000, 840000, -62400000, 887200000}}, 
  {{2102, b, 67264, 1}, {-54000, -5910000, 2577600000, -144012800000}}}, 
 {{{1109, a, 1109, 1}, {-504000, -170160000, -13802400000, -409512800000}}, 
  {{2218, a, 70976, 1}, {-826, 14724, -104799, 321556}}}, 
 {{{1197, a, 10773, 1}, {-410, 5316, -26047, 45716}}, {{1197, a, 410571, 1}, 
   {-670250, 343072500, -65325259375, 4388306862500}}, {{8379, a, 410571, 1}, 
   {-104000, -14910000, -372400000, -2312800000}}}, 
 {{{1462, a, 11696, 1}, {587750, -406017500, 80822490625, -4901824287500}}, 
  {{2924, a, 11696, 1}, {-120250, 8822500, 9615990625, 932688112500}}}, 
 {{{1532, a, 1532, 1}, {24750, -502500, 24365625, -18762500}}, 
  {{1532, a, 392192, 1}, {-205250, 6547500, 10087115625, -1046160012500}}, 
  {{124092, a, 372276, 1}, {-47250, 3217500, 141365625, 1533037500}}}, 
 {{{1757, a, 1757, 1}, {-14000, 1240000, -34400000, 1207200000}}, 
  {{12299, a, 602651, 1}, {-104000, -9590000, 607600000, -7487200000}}}, 
 {{{2016, c, 774144, 1}, {16000, -2040000, -38400000, -18547200000}}, 
  {{18144, d, 435456, 1}, {-24000, -360000, 213600000, 14227200000}}}, 
 {{{2301, a, 6903, 1}, {-15250, 1977500, -968634375, 87295637500}}, 
  {{62127, a, 62127, 1}, {-48, -496, -1152, -768}}}, 
 {{{3037, a, 3037, 1}, {72750, -8182500, 302365625, -3721962500}}, 
  {{6074, b, 194368, 1}, {-590250, 78352500, 53660740625, 4412311262500}}}, 
 {{{3161, a, 3161, 1}, {-15250, 897500, 41365625, 26237500}}, 
  {{256041, a, 768123, 1}, {-6000, 2170000, 39600000, -9669600000}}}, 
 {{{3434, a, 6868, 1}, {-10042, 635524, -15069023, 126933268}}, 
  {{6868, b, 439552, 1}, {-6440250, 10319102500, -6193808009375, 1320544367512500}}}, 
 {{{3434, b, 439552, 1}, {-64, 464, 2176, -13056}}, {{6868, a, 439552, 1}, 
   {-62250, 11892500, -946259375, 29696162500}}}, 
 {{{3456, b, 221184, 1}, {-72250, 2782500, 1331990625, -100889287500}}, 
  {{10368, a, 62208, 1}, {-24000, 2640000, 933600000, 57427200000}}}, 
 {{{3883, a, 469843, 1}, {-46000, 130000, 647600000, -20770400000}}, 
  {{28593, a, 85779, 1}, {-24000, 1610000, -86400000, 3772800000}}}, 
 {{{4698, b, 408726, 1}, {-57250, 1707500, 1889615625, 127767587500}}, 
  {{14094, a, 84564, 1}, {57750, -2507500, 804740625, -86233837500}}}, 
 {{{4968, a, 804816, 1}, {1044000, -217080000, 249609600000, -116481369600000}}, 
  {{14904, b, 89424, 1}, {-84000, 4040000, 1161600000, -39052800000}}}, 
 {{{5157, a, 5157, 1}, {-64, -112, 1152, 3328}}, {{15471, a, 417717, 1}, 
   {-66000, 1280000, 1401600000, -91490400000}}}, 
 {{{5496, a, 32976, 1}, {-76000, -3480000, 1273600000, 91750400000}}, 
  {{148392, a, 296784, 1}, {57750, -867500, 834740625, -24588037500}}}, 
 {{{6163, a, 6163, 1}, {-6000, -830000, -400000, 3530400000}}, 
  {{30815, a, 770375, 1}, {-5490345250, -257293893072500, -4521583643143384375, 
    -28252646059758163112500}}}, {{{7101, a, 7101, 1}, {26000, -360000, 153600000, -4072800000}}, 
  {{21303, a, 63909, 1}, {-54000, 3840000, 117600000, 2887200000}}}, 
 {{{7263, a, 21789, 1}, {6000, -190000, -50400000, 3532800000}}, 
  {{21789, b, 196101, 1}, {-410, -5116, -23103, -35884}}}, 
 {{{7680, a, 46080, 1}, {-36000, 680000, 305600000, -11110400000}}, 
  {{207360, a, 414720, 1}, {96, 64, 2304, 3072}}}, 
 {{{7680, b, 46080, 1}, {-36000, -680000, 305600000, 11110400000}}, 
  {{207360, b, 414720, 1}, {96, -64, 2304, -3072}}}, {{{8220, a, 295920, 1}, {-18, 92, -7, -836}}, 
  {{19180, a, 939820, 1}, {-21645250, 63688947500, -70273280884375, 27568762259987500}}}, 
 {{{8704, c, 591872, 1}, {-12702250, -25827807500, -18248673259375, -4113826556337500}}, 
  {{26112, c, 235008, 1}, {49750, -2127500, 68740625, -1003137500}}}, 
 {{{10154, a, 81232, 1}, {-136000, -6430000, 2575600000, -102789600000}}, 
  {{20308, a, 81232, 1}, {-50250, 7752500, -473759375, 12181262500}}}, 
 {{{11637, a, 11637, 1}, {-120250, -22197500, -1269009375, -23084987500}}, 
  {{34911, a, 34911, 1}, {-276000, -24230000, 20223600000, 3222650400000}}, 
  {{34911, a, 942597, 1}, {6000, 1690000, 9600000, 12067200000}}}, 
 {{{12150, a, 109350, 1}, {-1685250, 1016527500, -47638884375, -51247008112500}}, 
  {{24300, a, 291600, 1}, {-110250, -8527500, 1416740625, 975876862500}}}, 
 {{{12813, a, 115317, 1}, {-42250, 3132500, 422240625, -84513037500}}, 
  {{21355, a, 533875, 1}, {-37250, 6267500, -481884375, 16411787500}}}, 
 {{{13153, a, 644497, 1}, {26000, -1390000, 173600000, -627200000}}, 
  {{152199, b, 456597, 1}, {-2164000, 1965490000, -657884400000, 76998807200000}}}, 
 {{{13742, a, 54968, 1}, {-74000, -13360000, 1353600000, -28172800000}}, 
  {{34355, a, 858875, 1}, {437750, 3732500, 38216240625, 3133881962500}}}, 
 {{{15552, c, 746496, 1}, {16000, 2040000, -38400000, 18547200000}}, 
  {{62208, n, 746496, 1}, {-24000, 360000, 213600000, -14227200000}}}, 
 {{{17986, a, 17986, 1}, {-460250, -52777500, 55525490625, 11449983112500}}, 
  {{55545, a, 499905, 1}, {54, -204, -15, 420}}}, 
 {{{23808, f, 71424, 1}, {124000, 480000, 3453600000, -246950400000}}, 
  {{214272, a, 214272, 1}, {24000, 480000, 33600000, 6249600000}}}, 
 {{{49152, f, 294912, 1}, {-524000, 33640000, 60073600000, -13996172800000}}, 
  {{212992, c, 425984, 1}, {-14000, -1760000, 145600000, 13507200000}}}, 
 {{{61712, b, 987392, 1}, {-52250, -9557500, -784509375, -25437587500}}, 
  {{123424, b, 987392, 1}, {-21082, 1917380, -65025663, 780014356}}}};


(*################################################################*)

(* Section FOURtwo.  Examples involving Richelot isogenies. *)

richelotcheck[e_,f_,g_] :=
ABCDstarmat[
{-5*(7*e^2 - 2*f), -10*e*(3*e^2 - 2*f), 5*(32*e^4 - 39*e^2*f + g), 
 -4*e*(24*e^4 + 115*e^2*f - 5*g), -4*e*(80*e^4 + 7*e^2*f - g), 
 2*(40*e^4 - 9*e^2*f - g), -4*e*(5*e^2 + 2*f), 5*e^2 + 2*f}
]-
Module[{z=
5000*(125*e^4 + 20*f^2 - 4*g)^4*(25*e^2*f + g)^6},
{-5*(7*e^2 + 2*f), -10*e*(3*e^2 + 2*f), 5*(32*e^4 + 39*e^2*f + g), 
  -4*e*(24*e^4 - 115*e^2*f - 5*g)}*z^{2,3,4,5}]

(* The evaluation of ABCDstarmat[{-5*(7*e^2 - 2*f), -10*e*(3*e^2 - 2*f), 5*(32*e^4 - 39*e^2*f + g),
 -4*e*(24*e^4 + 115*e^2*f - 5*g), -4*e*(80*e^4 + 7*e^2*f - g),
 2*(40*e^4 - 9*e^2*f - g), -4*e*(5*e^2 + 2*f), 5*e^2 + 2*f}] referred to in 
the text was carried out with e=1 and took about half an hour; the evaluation 
for general then follows from homogeneity *)


(*################################################################*)

(* Section FOURthree.  Explicit families of modular abelian surfaces. *)

(* Many terms of A[a,b,c,d,s,t,u,v],...,D[a,b,c,d,s,t,u,v] vanish 
upon reduction modulo 3.  This allows symbolic computation of the 
discriminant, which is useful for checking the good reduction
condition of Corollary 3 *)

Delta3[{a_,b_,c_,d_,s_,t_,u_,v_}] := 
(2*a^3*b^2*c^2 + a^4*c^3 + a^2*c^4 + c^5 + a^3*b^3*d + 2*a^2*b*c^2*d + 
  2*b*c^3*d + 2*a*c^2*d^2 + 2*d^4)*
 (2*s^4 + 2*b*s*t^3 + c*t^4 + c*s*t^2*u + 2*c*s^2*u^2 + b^2*s*u^3 + 
   2*b*c*t*u^3 + 2*c^2*u^4 + 2*c*s^2*t*v + 2*a*c*t^3*v + a*c*s*t*u*v + 
   b*c*s*u^2*v + 2*c^2*t*u^2*v + a*b*c*u^3*v + a*c*s^2*v^2 + 2*b*c*s*t*v^2 + 
   c^2*t^2*v^2 + a^2*c*s*u*v^2 + c^2*s*u*v^2 + a*c^2*u^2*v^2 + a^3*b*s*v^3 + 
   b^3*s*v^3 + a*b*c*s*v^3 + 2*a^3*c*t*v^3 + 2*b^2*c*t*v^3 + a*c^2*t*v^3 + 
   2*b*c^2*u*v^3 + a^4*c*v^4 + a*b^2*c*v^4 + a^2*c^2*v^4 + c^3*v^4)^3*
 (2*s^12 + 2*b*s^9*t^3 + 2*a^3*b*s^3*t^9 + 2*b^3*s^3*t^9 + 2*a^3*b^2*t^12 + 
   2*b^4*t^12 + a^4*c*t^12 + b^2*s^9*u^3 + 2*a*c*s^9*u^3 + a^3*b^3*t^9*u^3 + 
   b^5*t^9*u^3 + 2*a*b^3*c*t^9*u^3 + a^3*c*d*t^9*u^3 + 2*a^3*b^4*s^3*u^9 + 
   b^6*s^3*u^9 + 2*a^3*c^3*s^3*u^9 + 2*b*d^3*s^3*u^9 + 2*a^3*b^5*t^3*u^9 + 
   b^7*t^3*u^9 + a^4*b^3*c*t^3*u^9 + 2*a^3*b*c^3*t^3*u^9 + 2*b^2*d^3*t^3*u^9 + 
   a*c*d^3*t^3*u^9 + a^3*b^6*u^12 + 2*b^8*u^12 + a*b^6*c*u^12 + 
   a^3*b^2*c^3*u^12 + 2*a^4*c^4*u^12 + a^3*b^3*c*d*u^12 + b^3*d^3*u^12 + 
   c*d^4*u^12 + a^3*b*s^9*v^3 + b^3*s^9*v^3 + c*d*s^9*v^3 + a^6*b^2*t^9*v^3 + 
   2*a^3*b^4*t^9*v^3 + b^6*t^9*v^3 + 2*a^7*c*t^9*v^3 + a^3*c^3*t^9*v^3 + 
   b^3*c*d*t^9*v^3 + a^6*b^5*u^9*v^3 + 2*b^9*u^9*v^3 + 2*a^7*b^3*c*u^9*v^3 + 
   a^6*b*c^3*u^9*v^3 + 2*a^3*b^3*c^3*u^9*v^3 + 2*b^6*c*d*u^9*v^3 + 
   a^3*c^4*d*u^9*v^3 + a^3*b^2*d^3*u^9*v^3 + b^4*d^3*u^9*v^3 + 
   2*a^4*c*d^3*u^9*v^3 + c^3*d^3*u^9*v^3 + a^12*b*s^3*v^9 + a^9*b^3*s^3*v^9 + 
   b^9*s^3*v^9 + 2*b*c^6*s^3*v^9 + b^4*d^3*s^3*v^9 + c^3*d^3*s^3*v^9 + 
   a^12*b^2*t^3*v^9 + a^9*b^4*t^3*v^9 + b^10*t^3*v^9 + 2*a^13*c*t^3*v^9 + 
   2*b^2*c^6*t^3*v^9 + a*c^7*t^3*v^9 + b^5*d^3*t^3*v^9 + 
   2*a*b^3*c*d^3*t^3*v^9 + b*c^3*d^3*t^3*v^9 + 2*a^12*b^3*u^3*v^9 + 
   2*a^9*b^5*u^3*v^9 + 2*b^11*u^3*v^9 + a^10*b^3*c*u^3*v^9 + a*b^9*c*u^3*v^9 + 
   b^3*c^6*u^3*v^9 + 2*a^12*c*d*u^3*v^9 + c^7*d*u^3*v^9 + 2*b^6*d^3*u^3*v^9 + 
   2*b^2*c^3*d^3*u^3*v^9 + a*c^4*d^3*u^3*v^9 + 2*b^3*c*d^4*u^3*v^9 + 
   2*a^15*b^2*v^12 + a^12*b^4*v^12 + 2*a^9*b^6*v^12 + 2*a^3*b^10*v^12 + 
   2*b^12*v^12 + a^16*c*v^12 + 2*a^12*c^3*v^12 + a^3*b^2*c^6*v^12 + 
   b^4*c^6*v^12 + 2*a^4*c^7*v^12 + c^9*v^12 + 2*a^9*b^3*c*d*v^12 + 
   2*b^9*c*d*v^12 + 2*a^3*b^5*d^3*v^12 + 2*b^7*d^3*v^12 + a^4*b^3*c*d^3*v^12 + 
   2*a^3*b*c^3*d^3*v^12 + b^3*c^3*d^3*v^12 + 2*c^4*d^4*v^12)^9

checkDelta3[{a_,b_,c_,d_,s_,t_,u_,v_}] := 
PolynomialMod[
Delta[PolynomialMod[ABCDnum[{a,b,c,d,s,t,u,v}],3]]-Delta3[{a,b,c,d,s,t,u,v}]
,3]

(* findisos finds the {s,t,u,v} that get one from {a,b,c,d} to {A,B,C,D} in the example: *)

modularexample := 
findisos[{12/5,12/5^2,292/5^3,-3672/5^5},{2^7/5,2^11*57/5^2,-2^12*503/5^3, 2^17*17943/5^5}]
(* the returned information includes the information {s,t,u,v} = {129/125, 11/25, 3/100, 1/20} 
appearing in the paper *)


(*################################################################*)

(* Section FOURfour.  Analogs for p=2 *)

traceless5[m_] := m - Tr[m]/5 IdentityMatrix[5]

newquintic[{a_,b_,c_,d_,s_,t_,u_,v_}] := 
Expand[
Module[{m},
m = {{0, 0, 0, 0, -d}, {1, 0, 0, 0, -c}, {0, 1, 0, 0, -b}, {0, 0, 1, 0, -a}, {0, 0, 0, 1, 0}};
-CharacteristicPolynomial[s m + t traceless5[m.m] + u traceless5[m.m.m] + v traceless5[m.m.m.m],x]]
]

(* The generic case newquintic[{a,b,c,d,s,t,u,v}] evaluates in about a twentieth of a second.  
In the resulting polynomial x^5+A x^3 + B x^2 + C x + D, the coefficients
A, B, C, and D respectively have 24, 86, 235, and 535 monomial terms, as asserted in the paper *)


(*################################################################*)
