/*
    Magma code related to the paper "Torsion subgroups of rational elliptic curves over the compositum of all cubic fields"
    by H. Daniels, A. Lozano-Robledo, F. Najman, and A.V. Sutherland
    Copyright (c) 2015 the authors.
    
    You are welcome to use and/or modify this code under GPL version 2 or any later version (but please be sure to cite the paper)

    Run VerifyAll() to verify all the computations cited in the paper (takes about 6 hours on a fast machine).
    The function Q3inftyTorsion(E) will compute E(Q(3^infty))_tor for any E/Q using the algorithm in described in Section 5.5.
    The function Q3inftyTorsionFast(E) is a faster version of the same thing, using the parameterizations in Table 1 (see Theorem 7.1).
*/


// Split Cartan -- diagonal matrices in GL(2,Z/nZ) -- E admits two distinct rational n-isogenies if and only if \im\rho_{E,n} is conjugate ot a subgroup of the split Cartan
function GL2SplitCartan(R)
    if Type(R) eq RngIntElt then R:=Integers(R); end if;
    M,pi:=MultiplicativeGroup(R);
    return sub<GL(2,R) | [[pi(g),0,0,1] : g in Generators(M)], [[1,0,0,pi(g)] : g in Generators(M)]>;
end function;

// Borel group -- upper triangular matrices in GL(2,Z/nZ) -- E admits a rational n-isogeny if and only if \im\rho_{E,n} is conjugate to a subgroup of  the Borel
function GL2Borel(R)
    if Type(R) eq RngIntElt then R:=Integers(R); end if;
    return sub<GL(2,R) | [1,1,0,1],GL2SplitCartan(R)>;
end function;

// Subgroup of Borel that fixes a point of order n
function GL2Borel1(R)
    if Type(R) eq RngIntElt then R:=Integers(R); end if;
    M,pi:=MultiplicativeGroup(R);    
    return sub<GL(2,R) | [1,1,0,1], [[1,0,0,pi(g)] : g in Generators(M)]>;
end function;

// Quick Monte-Carlo algorithm to check whether the Galois group of the splitting field of an integer poly
// has exponent not dividing 6 -- avoids wasting time computing the splitting field K of f
// If the algorithm returns true then K is not of generalized S3-type, if it returns false then it may or may not be.
function NotGenS3Poly(f)
    for i in [1..20] do
        while true do p:=RandomPrime(20); if p gt 100 then break; end if; end while;
        fp:=ChangeRing(f,GF(p));
        if fp eq 0 then continue; end if;
        if 6 mod LCM([Degree(g[1]):g in Factorization(fp)]) ne 0 then return true; end if;
    end for;
    return false;
end function;

// Fast algorithm to check whether a number field or group is of generalized S3-type via Lemma 3.2
function IsGenS3(G)
    if Type(G) eq FldRat then return true; end if; // as far as Magma is concerned, Q is not a number field
    if Type(G) eq FldNum then
        if NotGenS3Poly(DefiningPolynomial(G)) then return false; else return $$(GaloisGroup(G)); end if;
    end if;
    if Type(G) eq RngUPolElt then
        if NotGenS3Poly(G) then return false; else return $$(GaloisGroup(G)); end if;
    end if;
    return (6 mod Exponent(G)) eq 0 and IsAbelian(SylowSubgroup(G,2)) and  IsAbelian(SylowSubgroup(G,3)) and #{IsPrime(Index(G,H`subgroup)):H in MaximalSubgroups(G)} le 1;
end function;

// Checks whether H1 is conjugate to a subgroup of H2 in G
function IsConjugateToSubgroup(G,H1,H2)
    if not IsDivisibleBy(#H2,#H1) then return false; end if;
    n:=#H2 div #H1;
    return #[H:H in Subgroups(H2:IndexEqual:=n)|IsConjugate(G,H`subgroup,H1)] ne 0;
end function;

// Given two subgroups H1 and H2 of G, conjugate them to minimize the index of their intersection and return that intersection
function ConjugateIntersection(G,H1,H2)
        S:=[H1 meet H:H in Conjugates(G,H2)];
        I:=[<Index(G,S[i]),i>:i in [1..#S]];
        return S[Min(I)[2]];
end function;
    
// Given two subgroups H1 and H2 of G, conjugate them to minimize the size of the subgroup of G they generate
function ConjugateCompositum(G,H1,H2)
        S:=[sub<G|H1,H>:H in Conjugates(G,H2)];
        I:=[<#S[i],i>:i in [1..#S]];
        return S[Min(I)[2]];
end function;
    
// Returns the unique minimal normal subgroup H of G for which G/H is of generalized S3-type
function GenS3Divisor(G)
    if IsGenS3(G) then return sub<G|>; end if;
    S:=[H`subgroup:H in NormalSubgroups(G)];
    I:=Sort([<#S[i],i>:i in [1..#S]]);
    for i:=2 to #I do H:=S[I[i][2]]; if IsGenS3(quo<G|H>) then return H; end if; end for;
    assert false;    
end function;

// Compute the invariant factors of a finite Z/nZ-module of rank at most 2  (Magma apparently can't coerce Z/nZ-modules to abelian groups)
function ModuleInvariants(V)
    if Dimension(V) eq 0 then return []; end if;
    if Dimension(V) eq 1 then return [#V]; end if;
    assert Dimension(V) le 2;
    r1:=#sub<V|V.1>;  r2:=#sub<V|V.2>;
    return [GCD(r1,r2),LCM(r1,r2)];
end function;
   
// Given a subgroup of GL(2,Z/nZ), computes the invariants of the sub-module of Z/nZ x Z/nZ fixed by G (returns a list [], [a], or [a,b] with a|b|n)
function FixModule(H)
	V := Eigenspace(Identity(H),1);
	for h in Generators(H) do V:= V meet Eigenspace(Transpose(h),1); end for;	// take transpose to work with right eigenspaces
	return ModuleInvariants(V);
end function;

// Returns if the Z/nZ-module A contains a submodule isomorphic to B (where A and B are specified by invariants). 
function ModuleContains(A,B)
    i:=#A-#B;
    if i lt 0 then return false; end if;
    for j in [1..#B] do if not IsDivisibleBy(A[i+j],B[j]) then return false; end if; end for;
    return true;
end function;

// Given a subgroup G of GL(2,Z/nZ) corresponding to Gal(Q(E[n])/Q), returns [a,b] with a|b|n
// specifying the subgroup of Z/nZ x Z/nZ fixed by Gal(Q(E[n]) cap Q(3^infty))
function GenS3Torsion(G)
    return FixModule(GenS3Divisor(G));
end function;

// Returns true of a matrix group over a ring R contains elements of every possible determinant
function FullDeterminantMap(H)
    M,pi:=MultiplicativeGroup(BaseRing(H));
    return sub<M|[Inverse(pi)(Determinant(h)):h in Generators(H)]> eq M;
end function;

// Returns true if a give subgroup of GL(2,Z/nZ) contains an element corresponding to complex conjugation
// Here we use the stronger criterion from Remark 3.5 of  Sutherland's "Computing images of Galois representations attached to elliptic curves".
function ContainsComplexConjugation(H)
    return #[h:h in H|Determinant(h) eq -1 and Trace(h) eq 0 and ModuleContains(FixModule(sub<H|h>),[#BaseRing(H)])] gt 0;
end function;

// Checks that A is a possible set of invariants for  the torsion subgroup of an elliptic curve (i.e of the form Z/aZ x Z/bZ with a|b)
function ValidTorsionSubgroup(A)
    if #A gt 2 then return false; end if;
    if #A eq 0 then return true; end if;
    if #A eq 1 and A[1] gt 1 then return true; end if;
    return A[1] gt 1 and IsDivisibleBy(A[2],A[1]);
end function;

// Returns true if A contains a submodule isomorphic to B
function TorsionContains(A,B)
    if #B eq 0 or A eq B then return true; end if;
    assert ValidTorsionSubgroup(A) and ValidTorsionSubgroup(B);
    if #A lt #B then return false; end if;
    if not IsDivisibleBy(A[#A],B[#B]) then return false; end if;
    if #B eq 2 and not IsDivisibleBy(A[1],B[1]) then return false; end if;
    return true;
end function;

// Given a subgroup G of GL2(Z/nZ), returns the least divisor m for which G is the full inverse image of its reduction mod m
function GL2Level(G)
    if G eq GL(2,BaseRing(G)) then return 1; end if;
    return GCD([m:m in Divisors(#BaseRing(G)) | m gt 1 and Index(GL(2,Integers(m)),ChangeRing(G,Integers(m))) eq Index(GL(2,BaseRing(G)),G)]);
end function;

// given a torsion structure T with exponent n
// returns a list of the maximal subgroups G of GL(2,Z/nZ) for which the minimal normal subgroup H
// with G/H is generalized S3-type acts trivially on a subgroup of Z/nZ x Z/nZ isomorphic to T.
// These correspond to maximal mod-n Galois images of elliptic curves E/Q for which E(Q(3^infy))[n]
// is isomorphic to T.  The returned groups will each be reduced modulo m=GL2Level(G) (unless m=1)
function MaximalImagesWithGenS3Torsion(T)
    assert ValidTorsionSubgroup(T);
    assert #T eq 2;
    n := T[#T];
    A:=T;
    if n gt 2 and Valuation(n,2) eq 1 and #T eq 2 and Valuation(T[1],2) eq 1 and Valuation(T[2],2) eq 1 then
        print "invariants have 2-valuation 1, dividing by 2";
        n := n div 2; A:= [T[1] div 2, T[2] div 2];
        if A[1] eq 1 then A:=[A[2]]; end if;
    end if;
    s:=Cputime();
    assert IsDivisibleBy(n,A[#A]);
    for i:=1 to #A-1 do assert IsDivisibleBy(A[i+1],A[i]); end for;
    G:=GL(2,Integers(n));
    S0:=[H`subgroup:H in Subgroups(G)|FullDeterminantMap(H`subgroup) and ContainsComplexConjugation(H`subgroup)];
    S:=[H:H in S0|GenS3Torsion(H) eq A];
    t:=Cputime();
    printf "%o of %o candidate images G in GL(2,Z/%oZ) yielding torsion %o over Q(3^infty) found in %o secs, determining which are maximal...\n", #S, #S0, n, T, t-s;
    s:=Cputime();
    I:=Sort([<Index(G,S[i]),i>:i in [1..#S]]);
    indexes:={i[1]:i in I};
    U:=[<S[I[1][2]],I[1][1],AssociativeArray(indexes)>];
    for n in indexes do if IsDivisibleBy(n,U[1][2]) then U[1][3][n] := [K`subgroup: K in Subgroups(U[1][1]:IndexEqual:=n div U[1][2])]; end if; end for;
    for i:=2 to #I do
        H:=S[I[i][2]];
        Hn :=I[i][1];
        keep := true;
        for j:=1 to #U do
            if IsDivisibleBy(Hn,U[j][2]) then
                for K in U[j][3][Hn] do assert #H eq #K; if IsConjugate(G,H,K) then keep:= false; break; end if; end for;
                if not keep then break; end if;
            end if;
        end for;
        if keep then
            r:=<H,Hn,AssociativeArray(indexes)>;
            for n in indexes do if IsDivisibleBy(n,Hn) then r[3][n] := [K`subgroup:K in Subgroups(H:IndexEqual:=n div Hn)]; end if; end for;
            U:=Append(U,r);
        end if;
    end for;
    t:=Cputime();
    printf "%o maximal images G found for T = %o in %o secs\n", #U, T, t-s;
    S:=[r[1]:r in U];
    for i:= 1 to #S do m:=GL2Level(S[i]); if m gt 1 and m lt #BaseRing(S[i]) then S[i]:=ChangeRing(S[i],Integers(m)); end if; end for;
    return S;
end function;    

// given E/Q, determines if E(Q(3^infy))[p^n]=E[p^n]
// the cases p^n=7,9 are handled specially in order to keep the running time reasonable
function HasQ3inftyFullTorsion(E,p,n:UseSection5:=true,UseTheorem51:=true)
    assert IsPrime(p) and n ge 1;
    // Use Lemma 3.7 to rule out cases where Q(3^infty) does not contain any primitive p^nth roots of unity
    if p eq 5 or p^n gt 9 then return false; end if;
    // handle 7 specially by checking j(E) as given by Lemma 5.13
    if UseSection5 and p eq 7 then return jInvariant(E) eq 2268945/128; end if;
    // handle 9 specially by checking if j(E) is in the image of  one of the two possible j-maps permitted by Lemma 5.19 (j-maps taken from Sutherland-Zywina https://arxiv.org/abs/1605.03988)
    if UseSection5 and p^n eq 9 then
        F<t>:=FunctionField(Rationals());
        J3_7:=27*t^3*(8-t^3)^3/(t^3+1)^3;
        if #Roots(Numerator(J3_7-jInvariant(E))) gt 0 then return true; end if;
        J9_268:=432*t*(t^2-9)*(t^2+3)^3*(t^3-9*t+12)^3*(t^3+9*t^2+27*t+3)^3*(5*t^3-9*t^2-9*t-3)^3/((t^3-3*t^2-9*t+3)^9*(t^3+3*t^2-9*t-3)^3);
        if #Roots(Numerator(J9_268-jInvariant(E))) gt 0 then return true; end if;
        return false;
    end if;
    f,h:=HyperellipticPolynomials(E);
    assert h eq 0;
    g := DivisionPolynomial(E,p^n) div DivisionPolynomial(E,p^(n-1));
    for h in [a[1]:a in Factorization(g)] do
        if not IsGenS3(h) then return false; end if;
        K:=SplittingField(h);
        r:=Roots(ChangeRing(h,K))[1][1];
        R<u>:=PolynomialRing(K);
        L:=AbsoluteField(SplittingField(u^2-Evaluate(f,r)));
        if not IsGenS3(L) then return false; end if;
    end for;
    return true;
end function;

// given E/Q determines if E(Q(3^infy)) has point of order p^n
// the cases p^n=7,9 are handled specially in order to keep the running time reasonable
// by default we also use Theorem 5.1, but for verification this can be turned off.
function HasQ3inftyTorsionPoint(E,p,n:UseSection5:=true,UseTheorem51:=true)
    assert IsPrime(p) and n ge 1;
    // use Proposition 4.9 to rule out primes that we know are ruled out.
    if not p in [2,3,5,7,13] then return false; end if;
    // use Lemma 3.7, Thorem 4.4, and Lemma 4.6 to rule out powers of 7 and 13  (this just relies on the fact that no E/Q admits a 49-isogeny or a 169-isogeny)
    if p in [7,13] and n gt 1 then return false; end if;
    if UseTheorem51 then
        if p^n gt 42 or p^n in [25,27,32] then return false; end if;
    end if;
    // handle 7 specially by checking for special j(E) or if j(E) in in the image of the j-map for X0(7) (rational 7-isogeny), per Lemma 5.13
    if UseSection5 and p eq 7 and n eq 1 then
        F<t>:=FunctionField(Rationals());
        J7_76:=(t^2+5*t+1)^3*(t^2+13*t+49)/t;
        return jInvariant(E) eq 2268945/128 or #Roots(Numerator(J7_76-jInvariant(E))) gt 0;
    end if;
    // handle 9 specially by checking if j(E) is in the image of one of the two possible j-maps determined by Lemma 5.21 (j-maps taken from Sutherland-Zywina https://arxiv.org/abs/1605.03988)
    if UseSection5 and p eq 3 and n eq 2 then
        F<t>:=FunctionField(Rationals());
        J9_292:=(t+3)^3*(t^3+9*t^2+27*t+3)^3/(t*(t^2+9*t+27));
        if #Roots(Numerator(J9_292-jInvariant(E))) gt 0 then return true; end if;
        J9_299:=(t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3;
        if #Roots(Numerator(J9_299-jInvariant(E))) gt 0 then return true; end if;
        return false;
    end if;
    f,h:=HyperellipticPolynomials(E);
    assert h eq 0;
    g := DivisionPolynomial(E,p^n) div DivisionPolynomial(E,p^(n-1));
    for h in [a[1]:a in Factorization(g)] do
        if not IsGenS3(h) then continue; end if;
        K:=SplittingField(h);
        r:= Roots(ChangeRing(h,K))[1][1];
        R<u>:=PolynomialRing(K);
        L:=AbsoluteField(SplittingField(u^2-Evaluate(f,r)));
        if IsGenS3(L) then return true; end if;
    end for;
    return false;
end function;

// Returns the group invariants of the p-primary part of E(Q(3^infty)
function pPrimaryQ3inftyTorsion(E,p:UseSection5:=true,UseTheorem51:=true)
    E:=WeierstrassModel(IntegralModel(E));
    m:=0;
    while true do
        if not HasQ3inftyFullTorsion(E,p,m+1) then break; end if;
        m:=m+1;
    end while;
    n:=m;
    while true do
        if not HasQ3inftyTorsionPoint(E,p,n+1:UseSection5:=UseSection5,UseTheorem51:=UseTheorem51) then break; end if;
        n:=n+1;
    end while;
    return [p^m,p^n];
end function;

// Given an elliptic curve E (or a defining equation, Cremona referencee, or j-invariant != 1728) returns the invariants of the torsion subgroup of E(Q(3^infty))
// If UseTheorem51 is true (default) uses known upper bounds on the p-primary parts of E(Q(3^infty))_tor from Theorem 5.1
// If UseSection5 is true (default) uses Lemmas proved in Section 5 to speed up handling of 3-primary and 7-primary parts
function Q3inftyTorsion(E:UseSection5:=true,UseTheorem51:=true)
    if (Type(E) eq RngIntElt or Type(E) eq FldRatElt) then
        if E eq 1728 then print "For j-invariant 1728 you must specify a particular curve E/Q"; assert false; end if;
        E:=EllipticCurveFromjInvariant(E);
    end if;
    if Type(E) ne CrvEll then E:=EllipticCurve(E); end if;
    assert BaseRing(E) eq Rationals();
    T:=[1,1];  K:=Rationals();
    // Use Proposition 4.9
    for p in [2,3,5,7,13] do
        Tp:=pPrimaryQ3inftyTorsion(E,p:UseSection5:=UseSection5,UseTheorem51:=UseTheorem51);
        T[1]:=T[1]*Tp[1]; T[2]:=T[2]*Tp[2];
    end for;
    return T;
end function;

// maxN = 338 is sufficient to get every case except [14,14], to get them all, use maxN=2450
function Q3inftyTorsionExamples(maxN)
    examples:=AssociativeArray();
    D:=CremonaDatabase();
    for N in [11..maxN] do
        for E in EllipticCurves(D,N) do
           T:=Q3inftyTorsion(E);
            print CremonaReference(E),T;
            if IsDefined(examples,T) then
                examples[T]:=Append(examples[T],CremonaReference(E));
            else
                examples[T]:=[CremonaReference(E)];
            end if;
        end for;
    end for;
    return examples;
end function;
    
// Slow brute force implementation of Lemma 3.4 in http://www.math.cornell.edu/~zywina/papers/PossibleIndices/PossibleIndices.pdf
// For large n one will want to work with double cosets more intelligently, but this suffices for our purposes
function RationalCuspCount(H)
    GL2:=GL(2,BaseRing(H));
    U:=sub<GL2|[1,1,0,1],[-1,0,0,-1]>;
    M,pi:=MultiplicativeGroup(BaseRing(H));
    B:=sub<GL2|[[pi(g),0,0,1]:g in Generators(M)]>;
    // compute right cosets H\GL2
    R:={{h:h in H}};
    for g in GL2 do
        if #{C:C in R|g in C} eq 0 then R := R join {{h*g:h in H}}; end if;
    end for;
    // compute double costs H\GL2/U
    LR:={{s*u:u in U, s in S}:S in R};
     // return # of double cosets fixed by the right-action of B, along with the total number
    return #[S:S in LR|{s*b:b in B, s in S} eq S];
end function;
    
// Groups taken from Cummins-Pauli tables at http://www.uncg.edu/mat/faculty/pauli/congruence/congruence.html
CPGroups:=AssociativeArray();
CPGroups["1A0"]:=sub<SL(2,Integers(2))|[1,1,0,1],[0,1,1,0]>;    // We treat 1A0 as a subgroup of SL2(Z/2Z) (no matrix groups over the ring with one element)
CPGroups["2C0"]:=sub<SL(2,Integers(2))|>;
CPGroups["3B0"]:=sub<SL(2,Integers(3))|[2,1,2,0],[0,2,1,2],[1,2,1,0],[2,0,0,2]>;
CPGroups["3D0"]:=sub<SL(2,Integers(3))|[2,0,0,2]>;
CPGroups["4B0"]:=sub<SL(2,Integers(4))|[0,1,3,2],[2,1,3,0],[1,2,2,1],[3,2,2,3]>;
CPGroups["4D0"]:=sub<SL(2,Integers(4))|[3,0,0,3],[1,1,1,2]>;
CPGroups["4E0"]:=sub<SL(2,Integers(4))|[1,2,2,1],[3,2,2,3]>;
CPGroups["4G0"]:=sub<SL(2,Integers(4))|[3,0,0,3]>;
CPGroups["5D0"]:=sub<SL(2,Integers(5))|[2,3,3,0],[4,0,0,4]>;
CPGroups["6I0"]:=sub<SL(2,Integers(6))|[5,4,2,3],[3,4,2,1],[3,2,4,5]>;
CPGroups["7B0"]:=sub<SL(2,Integers(7))|[3,2,0,5],[6,0,0,6],[5,3,4,4]>;
CPGroups["7B1"]:=sub<SL(2,Integers(7))|[4,0,2,2],[6,0,0,6]>;
CPGroups["8C0"]:=sub<SL(2,Integers(8))|[5,4,4,5],[1,4,4,1],[4,5,3,2],[7,2,6,3],[3,0,0,3]>;
CPGroups["8G0"]:=sub<SL(2,Integers(8))|[1,2,0,1],[5,4,0,5],[1,4,0,1],[1,6,0,1],[3,2,0,3]>;
CPGroups["8O0"]:=sub<SL(2,Integers(8))|[7,4,0,7],[1,4,0,1],[5,6,0,5]>;
CPGroups["9B0"]:=sub<SL(2,Integers(9))|[4,0,6,7],[7,0,6,4],[8,0,0,8],[7,0,8,4]>;
CPGroups["9C0"]:=sub<SL(2,Integers(9))|[4,0,6,7],[7,0,6,4],[1,3,2,7],[8,0,0,8]>;
CPGroups["9H0"]:=sub<SL(2,Integers(9))|[4,3,6,7],[4,0,0,7],[8,0,0,8]>;
CPGroups["9I0"]:=sub<SL(2,Integers(9))|[1,0,6,1],[8,0,0,8],[7,0,5,4]>;
CPGroups["9J0"]:=sub<SL(2,Integers(9))|[1,0,6,1],[1,3,2,7],[8,0,0,8]>;
CPGroups["9A1"]:=sub<SL(2,Integers(9))|[4,0,6,7],[7,0,6,4],[8,0,0,8],[4,6,5,1]>;
CPGroups["9C1"]:=sub<SL(2,Integers(9))|[1,0,6,1],[8,0,0,8],[4,0,3,7]>;
CPGroups["9H1"]:=sub<SL(2,Integers(9))|[1,0,3,1],[8,0,0,8]>;
CPGroups["12E0"]:=sub<SL(2,Integers(12))|[7,0,6,7],[1,0,6,1],[9,4,8,1],[9,8,1,5],[11,4,8,3],[9,8,4,5],[11,4,5,3]>;
CPGroups["12I1"]:=sub<SL(2,Integers(12))|[10,3,3,1],[7,0,0,7],[5,0,0,5],[1,0,8,1],[11,0,8,11],[5,0,4,5]>;
CPGroups["12P1"]:=sub<SL(2,Integers(12))|[5,4,8,9],[7,0,0,7],[1,8,4,9],[9,8,4,5],[9,8,10,5]>;
CPGroups["13B0"]:=sub<SL(2,Integers(13))|[4,0,11,10],[0,3,4,2],[5,7,3,7]>;
CPGroups["14E2"]:=sub<SL(2,Integers(14))|[3,0,6,5],[3,10,8,13],[11,8,12,5],[13,10,2,7]>;
CPGroups["15G1"]:=sub<SL(2,Integers(15))|[7,3,3,10],[4,0,10,4],[1,0,5,1],[11,0,5,11]>;
CPGroups["16H0"]:=sub<SL(2,Integers(16))|[13,12,4,5],[2,1,15,0],[13,6,10,1],[9,8,8,9],[13,4,12,5],[15,6,10,3],[1,8,8,1]>;
CPGroups["16G1"]:=sub<SL(2,Integers(16))|[13,12,4,5],[13,6,10,1],[0,1,15,6],[9,8,8,9],[15,6,10,3],[1,8,8,1],[2,3,13,4]>;
CPGroups["16M1"]:=sub<SL(2,Integers(16))|[1,12,0,1],[13,10,8,5],[13,2,8,5],[15,4,0,15],[1,8,0,1],[9,12,0,9]>;
CPGroups["21B1"]:=sub<SL(2,Integers(21))|[13,0,0,13],[1,0,14,1],[8,0,0,8],[1,0,7,1],[1,0,18,1],[20,0,20,20],[16,0,0,4]>;
CPGroups["27A0"]:=sub<SL(2,Integers(27))|[26,0,0,26],[10,0,24,19],[10,0,0,19],[1,0,18,1],[13,9,1,7],[16,18,25,13]>;
CPGroups["27A2"]:=sub<SL(2,Integers(27))|[26,0,0,26],[10,0,24,19],[22,0,0,16],[10,0,0,19],[1,0,18,1],[13,9,1,7]>;
CPGroups["27B4"]:=sub<SL(2,Integers(27))|[10,0,15,19],[1,0,15,1],[10,0,0,19],[1,18,2,10],[1,0,18,1],[17,0,0,8]>;
CPGroups["27A4"]:=sub<SL(2,Integers(27))|[10,0,15,19],[1,0,15,1],[10,0,0,19],[7,9,19,13],[1,0,18,1],[17,0,0,8]>;
CPGroups["27C1"]:=sub<SL(2,Integers(27))|[26,0,0,26],[10,0,0,19],[19,0,24,10],[25,9,20,4],[1,0,18,1]>;
CPGroups["28A2"]:=sub<SL(2,Integers(28))|[13,24,16,21],[15,0,0,15],[15,14,7,1],[25,8,12,5],[1,7,7,22],[13,12,12,9]>;
CPGroups["28E3"]:=sub<SL(2,Integers(28))|[13,0,0,13],[17,24,8,13],[13,24,16,21],[15,0,0,15],[1,7,7,22]>;
CPGroups["28E5"]:=sub<SL(2,Integers(28))|[19,8,18,15],[1,0,14,1],[25,16,4,25],[5,8,18,1],[25,8,12,5],[5,8,4,1]>;

// Referneced groups from Rouse and Zureick-Brown "Images of the 2-adic Galois representation attached to elliptic curves", data from http://users.wfu.edu/rouseja/2adic/
RZBGroups:=AssociativeArray();
RZBGroups[["H20"]]:=sub<GL(2,Integers(4))|[3,3,0,1],[3,0,0,3],[1,1,3,0]>;
RZBGroups[["H335"]]:=sub<GL(2,Integers(16))|[9,0,0,1],[1,8,0,1],[15,0,0,15],[15,10,0,15],[11,9,8,13],[9,2,0,9],[7,8,0,7],[3,15,8,3],[15,4,0,15],[15,4,0,3],[1,12,0,1]>;
RZBGroups[["H478"]]:=sub<GL(2,Integers(16))|[1,12,0,1],[15,4,0,3],[13,5,8,15],[7,10,0,11],[15,0,0,15],[9,12,0,1],[7,8,0,7],[1,8,0,1]>;

// Maps to X(1) taken from Sutherland-Zywina "Modular curves of prime-power level with infinitely many rational points" https://arxiv.org/abs/1605.03988.
F<t> := FunctionField(Rationals());
SZGroups:=AssociativeArray();
SZGroups[[2,0]]:=<"2C0",sub<GL(2,Integers(2))|>,(t^2+192)^3/(t^2 - 64)^2>;
SZGroups[[3,9]]:=<"3B0",sub<GL(2,Integers(3))|[2,1,2,0],[1,0,1,2],[2,0,0,2]>,(t+3)^3*(t+27)/t>;
SZGroups[[3,7]]:=<"3D0",sub<GL(2,Integers(3))|[2,0,0,1],[2,0,0,2]>, 27*t^3*(8-t^3)^3/(t^3+1)^3>;
SZGroups[[4,55]]:=<"4A0",sub<GL(2,Integers(4))|[3,1,0,1],[3,0,0,3],[0,3,1,3],[3,2,2,1]>,-4*t^4 - 32*t^3>;
SZGroups[[4,46]]:=<"4B0",sub<GL(2,Integers(4))|[3,0,0,1],[1,1,0,1],[3,0,0,3],[1,2,0,1]>,(t^2+16*t+16)^3/(t*(t+16))>;
SZGroups[[4,42]]:=<"4D0",sub<GL(2,Integers(4))|[3,1,0,1],[3,0,0,3],[0,3,1,3]>,-4*(t^2 + 2*t - 2)^3*(t^2 + 10*t - 2)/t^4>;
SZGroups[[4,32]]:=<"4E0",sub<GL(2,Integers(4))|[3,1,0,1],[3,0,0,3],[1,2,0,1]>,-16*(t^4-14*t^2+1)^3/(t^2*(t^2+1)^4)>;
SZGroups[[4,30]]:=<"4E0",sub<GL(2,Integers(4))|[3,0,0,1],[3,0,0,3],[1,2,0,1]>,16*(t^4+4*t^3+20*t^2+32*t+16)^3/(t^4*(t+1)^2*(t+2)^4)>;
SZGroups[[4,15]]:=<"4G0",sub<GL(2,Integers(4))|[3,0,0,1],[3,0,0,3]>,(t^8+224*t^4+256)^3/(t^4*(t^4-16)^4)>;
SZGroups[[5,38]]:=<"5D0",sub<GL(2,Integers(5))|[2,1,2,3],[4,1,4,4],[0,3,3,2],[4,0,0,4]>,
(t^4-12*t^3+14*t^2+12*t+1)^3/(t^5*(t^2-11*t-1))>;
SZGroups[[7,76]]:=<"7B0",sub<GL(2,Integers(7))|[5,6,3,2],[2,5,4,0],[6,4,2,4],[6,0,0,6],[2,0,0,2]>,
(t^2+5*t+1)^3*(t^2+13*t+49)/t>;
SZGroups[[8,2142]]:=<"8C0",sub<GL(2,Integers(8))|[1,1,0,1],[1,2,0,1],[3,1,0,5],[5,0,0,5],[1,4,0,1],[3,0,0,1],[3,0,0,3]>,
(t^4-16*t^2+16)^3/(t^2*(t^2-16))>;
SZGroups[[8,1942]]:=<"8G0",sub<GL(2,Integers(8))|[1,2,0,1],[3,1,0,5],[5,0,0,5],[3,1,0,1],[1,4,0,1],[3,0,0,3]>,
-4*(t^8-60*t^6+134*t^4-60*t^2+1)^3/(t^2*(t^2-1)^2*(t^2+1)^8)>;
SZGroups[[8,1545]]:=<"8O0",sub<GL(2,Integers(8))|[1,0,0,5],[1,2,0,1],[7,0,0,7],[1,0,0,3],[1,4,0,1]>,
(t^16-8*t^14+12*t^12+8*t^10+230*t^8+8*t^6+12*t^4-8*t^2+1)^3/(t^8*(t^2-1)^8*(t^2+1)^4*(t^4- 6*t^2 + 1)^2)>;
SZGroups[[9,292]]:=<"9B0",sub<GL(2,Integers(9))|[1,1,0,1],[1,3,0,1],[4,0,0,1],[8,0,0,1],[8,0,0,8],[4,0,0,4]>,
(t+3)^3*(t^3+9*t^2+27*t+3)^3/(t*(t^2+9*t+27))>;
SZGroups[[9,299]]:=<"9C0",sub<GL(2,Integers(9))|[1,3,0,1],[4,0,0,1],[8,0,0,1],[4,5,3,1],[8,0,0,8],[4,0,0,4]>,
(t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3>;
SZGroups[[9,259]]:=<"9H0",sub<GL(2,Integers(9))|[1,3,0,1],[4,0,0,1],[7,6,3,4],[8,7,0,1],[8,0,0,8]>,
1/8*(t^3-3*t^2-9*t+3)^3*(t^3+9*t^2-9*t-9)^3*(t^6-18*t^5+171*t^4+180*t^3-297*t^2-162*t+189)^3/((t^2-1)^3*(t^2+3)^9*(t^3-9*t^2-9*t+9)^3)>;
SZGroups[[9,245]]:=<"9I0",sub<GL(2,Integers(9))|[1,1,0,1],[1,3,0,1],[8,0,0,1],[8,0,0,8],[1,0,0,4]>,
(t^3-3*t^2+1)^3*(t^9-9*t^8+27*t^7-48*t^6+54*t^5-45*t^4+27*t^3-9*t^2+1)^3/(t^9*(t-1)^9*(t^2-t+1)^3*(t^3-6*t^2+3*t+1))>;
SZGroups[[9,275]]:=<"9J0",sub<GL(2,Integers(9))|[1,3,0,1],[8,0,0,1],[4,5,3,1],[8,0,0,8],[1,0,0,4]>,
(t^2-t+1)^3*(t^3+3*t^2-6*t+1)*(t^9-9*t^7+6*t^6+18*t^5-9*t^4-27*t^3+27*t^2-9*t+1)^3/(t^9*(t-1)^9*(t^3-3*t+1)^3)>;
SZGroups[[9,268]]:=<"9J0",sub<GL(2,Integers(9))|[1,3,0,1],[4,5,3,1],[8,7,0,1],[8,0,0,8],[4,0,0,4]>,
432*t*(t^2-9)*(t^2+3)^3*(t^3-9*t+12)^3*(t^3+9*t^2+27*t+3)^3*(5*t^3-9*t^2-9*t-3)^3/((t^3-3*t^2-9*t+3)^9*(t^3+3*t^2-9*t-3)^3)>;
SZGroups[[13,204]]:=<"13B0",sub<GL(2,Integers(13))|[3,11,4,10],[1,0,5,5],[9,0,0,9],[12,0,0,12],[0,7,11,2],[2,7,12,10]>,
(t^4-t^3+5*t^2+t+1)*(t^8-5*t^7+7*t^6-5*t^5+5*t^3+7*t^2+5*t+1)^3/(t^13*(t^2-3*t-1))>;
SZGroups[[16,25902]]:=<"16H0",sub<GL(2,Integers(16))|[1,1,0,1],[9,0,0,1],[1,2,0,1],[9,0,0,3],[7,0,0,7],[9,0,0,9],[1,4,0,1],[1,8,0,1],[9,11,0,5]>,
(t^16-8*t^14+12*t^12+8*t^10-10*t^8+8*t^6+12*t^4-8*t^2+1)^3/((t^4-6*t^2+1)*(t^2+1)^2*(t^2-1)^4*t^16)>;
SZGroups[[27,2920]]:=<"27A0",sub<GL(2,Integers(27))|[1,1,0,1],[1,0,3,26],[26,0,0,26],[10,0,0,10],[1,3,0,1],[16,0,18,4],[4,0,0,4],[10,0,0,1],[1,9,0,1]>,
(t^3+3)^3*(t^9+9*t^6+27*t^3+3)^3/(t^3*(t^6+9*t^3+27))>;

// Separate verification procedures for each Example/Lemma/Remark/Proposition/Theorem appear below -- use VerifyAll() to run all of them
// Note that results that rely purely on citations to the literature (or pure thought) are not verified here, only those that reference a computation or check

procedure VerifyExample19()
    print "This will take several minutes...";
    A:=Q3inftyTorsionExamples(338);
    E:=MinimalTwist(EllipticCurveFromjInvariant(2268945/128));
    A[Q3inftyTorsion(E)]  := [CremonaReference(E)];
    assert #A eq 20;
    assert A[[2,2]][1] eq "11a2";
    assert A[[2,4]][1] eq "17a3";
    assert A[[2,8]][1] eq "15a5";
    assert A[[2,10]][1] eq "11a1";
    assert A[[2,14]][1] eq "26b1";
    assert A[[2,16]][1] eq "210e1";
    assert A[[2,26]][1] eq "147b1";
    assert A[[4,4]][1] eq "17a1";
    assert A[[4,8]][1] eq "15a2";
    assert A[[4,16]][1] eq "210e2";
    assert A[[4,28]][1] eq "338a1";
    assert A[[6,6]][1] eq "20a1";
    assert A[[6,12]][1] eq "30a1";
    assert A[[6,18]][1] eq "14a3";
    assert A[[6,30]][1] eq "50a3";
    assert A[[6,42]][1] eq "162b1";
    assert A[[8,8]][1] eq "15a1";
    assert A[[12,12]][1] eq "30a2";
    assert A[[14,14]][1] eq "2450a1";
    assert A[[18,18]][1] eq "14a1";
end procedure;

// Note: for the results in section 5, calls to Q3inftyTorsion should set UseSection5 to false

procedure VerifyLemma510()
    for G in [H`subgroup:H in Subgroups(GL(2,Integers(8)))| FullDeterminantMap(H`subgroup)] do
        if IsGenS3(G) then assert Exponent(G) eq 2 and IsAbelian(G); end if;
    end for;
end procedure;

procedure VerifyCorollary515()
    assert jInvariant(EllipticCurve("2450a1")) eq 2268945/128;
    assert Q3inftyTorsion(EllipticCurve("2450a1"):UseSection5:=false) eq [14,14];
end procedure;
    
procedure VerifyLemma516()
    G3:=GL(2,Integers(3));
    B:=sub<G3|[1,0,0,2],[2,0,0,1],[1,1,0,1]>;
    for G in [G`subgroup : G in Subgroups(G3)] do
        assert IsConjugateToSubgroup(G3,G,B) eq IsGenS3(G); 
    end for;
end procedure;

procedure VerifyLemma517()
    print "This will take a few minutes...";
    S:=MaximalImagesWithGenS3Torsion([9,27]);
    assert #S eq 1;
    G9:=GL(2,Integers(9));
    G3,pi:=ChangeRing(G9,Integers(3));
    S9:=SL(2,Integers(9));
    S3:=pi(S9);
    H:=sub<G9|[1,3,0,1],[8,0,0,1],[1,0,0,2]>;
    assert IsConjugate(G9,S[1],H);
    assert IsConjugate(S9,H meet S9, CPGroups["9H1"]);
    H1:=SZGroups[[9,245]][2];  H2:=SZGroups[[9,275]][2];
    h1:=SZGroups[[9,245]][3]; h2:=SZGroups[[9,275]][3];
    assert IsConjugate(S9,S9 meet H1, CPGroups["9I0"]);
    assert IsConjugate(S9,S9 meet H2, CPGroups["9J0"]);
    B:=SZGroups[[3,9]][2]; b:=SZGroups[[3,9]][3];
    assert IsConjugate(G3,B,GL2Borel(3));
    assert IsConjugate(G9,H,ConjugateIntersection(G9,H1,H2));
    assert IsConjugate(G9,sub<G9|Kernel(pi),Inverse(pi)(B)>,ConjugateCompositum(G9,H1,H2));
    // maps from X_H1, X_H2 to X_0(3), from SZ
    F<t>:=FunctionField(Rationals());
    f1:=(t^2-t+1)^3*(t^3-6*t^2+3*t+1)/(t^3*(t-1)^3);
    f2:=(t^3-3*t+1)^3/(t^3*(t-1)^3);
    assert Evaluate(b,f1) eq h1 and Evaluate(b,f2) eq h2;
    // compute a model for X_H as the fiber product over X_0(3) of the maps from X_H1, X_H2, and X_H3
    R<x,y>:=PolynomialRing(Rationals(),2);
    g:=Numerator(Evaluate(f1,x)-Evaluate(f2,y));
    // the polynomial g has 3 factors, each of which defines a genus 1 curve isomorphic to 27a3
    // the issue here is that X_H1 and X_H2 don't exactly cut out X_H.  There is an X_H3 we could add
    // (using the group H3=[9,259]) to make this true, but then would would have to work in P^3 and it gets too time consuming
    for f in [f[1]:f in Factorization(g)] do
        C:=ProjectiveClosure(Curve(AffineSpace(R),f));
        assert Genus(C) eq 1;
        P:=[P:P in RationalPoints(C:Bound:=100)|not P in SingularPoints(C)][1];
        E:=EllipticCurve(C,P);
        assert CremonaReference(E) eq "27a3";
    end for;
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 3;
    assert RationalCuspCount(H) eq 3;
end procedure;
    
procedure VerifyLemma518()
    print "This will take a few minutes...";
    S:=MaximalImagesWithGenS3Torsion([3,27]);
    assert #S eq 3;
    G27:=GL(2,Integers(27));
    G9,pi:=ChangeRing(G27,Integers(9));
    S27:=SL(2,Integers(27));
    S9:=SL(2,Integers(9));
    H1:=[H:H in S|IsConjugate(S27,H meet S27,CPGroups["27C1"])][1];
    H2:=[H:H in S|IsConjugate(S27,H meet S27,CPGroups["27B4"])][1];
    H3:=[H:H in S|IsConjugate(S27,H meet S27,CPGroups["27A4"])][1];
    assert IsConjugate(G27,H3,sub<G27|[1,2,9,1],[1,0,0,2],[8,0,0,1]>);
    // check that H1 is conjugate to a subgroup of the Borel
    assert IsConjugateToSubgroup(G27,H1,GL2Borel(27));
    // check that none of the curves in isogeny class 27a have a point of order 27 over Q(3^infty)
    for E in EllipticCurves(CremonaDatabase(),27) do
        assert jInvariant(E) ne 1728;
        f:=DivisionPolynomial(E,27) div DivisionPolynomial(E,9);
        for g in [a[1]:a in Factorization(f)] do
            assert not IsGenS3(g);
        end for;
    end for;
    H4:=sub<G27|[1,1,9,1],[1,0,0,2],[2,0,0,1]>;
    assert IsConjugate(S27,H4 meet S27, CPGroups["27A2"]);
    assert IsConjugateToSubgroup(G27,H2,H4);
    // show that the modular curve X_{H4} which is a quotient of X_{H2} has no non-cuspidal rational points
    assert RationalCuspCount(H4) eq 2;
    R<x>:=PolynomialRing(Rationals());
    C:=HyperellipticCurve(x^6-18*x^3-27);
    assert #Chabauty0(Jacobian(C)) eq 2;
    assert RationalCuspCount(H4) eq 2;
    // now construct our model for X_{H_3}
    // the groups and maps below are from [SZ] and only included for reference
    H5:=SZGroups[[27,2920]][2]; H6:=SZGroups[[9,245]][2];
    h5:=SZGroups[[27,2920]][3]; h6:=SZGroups[[9,245]][3];  
    assert IsConjugate(S27,H5 meet S27,CPGroups["27A0"]);
    assert IsConjugate(S9,H6 meet S9,CPGroups["9I0"]);
    B:=SZGroups[[9,292]][2]; b:=SZGroups[[9,292]][3];
    // this shows that X_H3 can be constructed as a fiber product over the maps from X_H5 and X_H6 to X_B=X_0(9)
    assert IsConjugate(G27,H3,ConjugateIntersection(G27,H5,sub<G27|Kernel(pi),Inverse(pi)(H6)>));
    assert IsConjugate(G9,GL2Borel(9),B);
    assert IsConjugate(G27,sub<G27|Kernel(pi),Inverse(pi)(B)>,ConjugateCompositum(G27,H5,sub<G27|Kernel(pi),Inverse(pi)(H6)>));
    // maps from X_H5 and X_H6 to X_0(9)
    F<t>:=FunctionField(Rationals());
    f5 := t^3;
    f6 := (t^3-6*t^2+3*t+1)/(t^2-t);
    assert Evaluate(b,f5) eq h5 and Evaluate(b,f6) eq h6;
    R<x,y>:=PolynomialRing(Rationals(),2);
    f:=Numerator(Evaluate(f5,x)-Evaluate(f6,y));
    assert IsIrreducible(f) and f eq x^3*y^2-x^3*y-y^3+6*y^2-3*y-1;
    C:=ProjectiveClosure(Curve(AffineSpace(R),f));
    assert Genus(C) eq 4;
    assert #SingularPoints(C) eq 2;
    assert RationalPoints(C:Bound:=1000) eq SingularPoints(C);
    // see separate Magma script Q3torGenus4Case.m by Jackson Morrow for the analysis of the rational points on C
end procedure;
        
procedure VerifyLemma519()
    print "This will take a minute...";
    G9:=GL(2,Integers(9));
    S9:=SL(2,Integers(9));
    G3,pi:=ChangeRing(G9,Integers(3));
    S3:=pi(S9);
    C:=GL2SplitCartan(Integers(3));
    H:=sub<G9|[1,2,3,1],[1,3,0,1],[1,0,0,8],[2,0,0,2]>;
    assert IsConjugate(G3,C,SZGroups[[3,7]][2]);
    assert IsConjugate(G9,H,SZGroups[[9,268]][2]);
    assert IsConjugate(S3,C meet S3,CPGroups["3D0"]);
    assert IsConjugate(S9,H meet S9,CPGroups["9J0"]);
    F<t>:=FunctionField(Rationals());
    assert 27*t^3*(8-t^3)^3/(t^3+1)^3 eq SZGroups[[3,7]][3];
    assert 432*t*(t^2-9)*(t^2+3)^3*(t^3-9*t+12)^3*(t^3+9*t^2+27*t+3)^3*(5*t^3-9*t^2-9*t-3)^3/((t^3-3*t^2-9*t+3)^9*(t^3+3*t^2-9*t-3)^3) eq SZGroups[[9,268]][3];
    C9:=sub<G9|Kernel(pi),Inverse(pi)(C)>;
    assert IsGenS3(C9) and IsGenS3(H);
    for G in [G`subgroup:G in Subgroups(G9)|IsGenS3(G`subgroup)] do
        assert IsConjugateToSubgroup(G9,G,C9) or IsConjugateToSubgroup(G9,G,H);
    end for;
end procedure;

procedure VerifyExample520()
    assert pPrimaryQ3inftyTorsion(EllipticCurve("27a3"),3:UseSection5:=false) eq [9,9];
    assert pPrimaryQ3inftyTorsion(EllipticCurve("17100g2"),3:UseSection5:=false) eq [9,9];
end procedure;

procedure VerifyLemma521()
    print "This will take a minute...";
    F<t>:=FunctionField(Rationals());
    G9:=GL(2,Integers(9));
    S9:=SL(2,Integers(9));
    G3,pi:=ChangeRing(G9,Integers(3));
    H1:=sub<G9|[1,1,0,1],[2,0,0,1],[2,0,0,2]>;
    H2:=sub<G9|[1,2,3,1],[2,0,0,1],[2,0,0,2]>;
    H3:=sub<G9|[1,1,3,1],[2,0,0,1],[2,0,0,2]>;
    assert IsConjugate(G9,H1,SZGroups[[9,292]][2]);
    assert (t+3)^3*(t^3+9*t^2+27*t+3)^3/(t*(t^2+9*t+27)) eq SZGroups[[9,292]][3];
    assert (t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3 eq SZGroups[[9,299]][3];
    assert IsConjugate(G9,H2,SZGroups[[9,299]][2]);
    assert IsConjugate(S9,H3 meet S9,CPGroups["9A1"]);
    T9:=sub<G9|[1,1,0,1],[2,0,0,1]>;
    assert not IsGenS3(H1) and not IsGenS3(H2) and not IsGenS3(H3);
    N1:=GenS3Divisor(H1);
    assert IsGenS3(quo<H1|N1>) and #pi(N1) eq 1 and IsConjugateToSubgroup(G9,N1,T9);
    N2:=GenS3Divisor(H2);
    assert IsGenS3(quo<H2|N2>) and #pi(N2) eq 1 and IsConjugateToSubgroup(G9,N2,T9);
    N3:=GenS3Divisor(H3);
    assert IsGenS3(quo<H3|N3>) and #pi(N3) eq 1 and IsConjugateToSubgroup(G9,N3,T9);
    S:=[G`subgroup:G in Subgroups(G9)|not IsGenS3(G`subgroup) and #pi(N) eq 1 and IsConjugateToSubgroup(G9,N,T9) where N:=GenS3Divisor(G`subgroup)];
    for G in S do
        assert IsConjugateToSubgroup(G9,G,H1) or IsConjugateToSubgroup(G9,G,H2) or IsConjugateToSubgroup(G9,G,H3);
    end for;
    E:=EllipticCurve("27a3");   // model for X_H3 obtained by RZB
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 3;
    assert RationalCuspCount(H3) eq 2;
    assert pPrimaryQ3inftyTorsion(EllipticCurveFromjInvariant(0),3:UseSection5:=false) eq [9,9];
end procedure;
        
procedure VerifyLemma61()
    G4:=GL(2,Integers(4));
    G2,pi:=ChangeRing(G4,Integers(2));
    S:=[H:H in [G`subgroup:G in Subgroups(G4)] | FullDeterminantMap(H) and ContainsComplexConjugation(H) and #pi(H) ge 3 and IsGenS3(H)];
    H:=sub<G4|[0,3,1,3],[3,1,0,1],[3,0,0,3]>; // this is SZGroups([[4,42]]), but we don't cite SZ in Lemma 6.1 so we won't verify it here
    for G in S do assert IsConjugateToSubgroup(G4,G,H); end for;
    F<t>:=FunctionField(Rationals());
    jt:=-4*(t^2-3)^3*(t^2-8*t-11)/(t+1)^4;
    assert #Roots(Numerator(jt-0)) eq 0 and #Roots(Numerator(jt-1728)) eq 0;
    E:=EllipticCurveFromjInvariant(jt);
    assert IsSquare((jt-1728)/Discriminant(E)) and IsSquare(-Discriminant(E));
end procedure;

procedure VerifyProposition62()
    F<t>:=FunctionField(Rationals());
    jt:=(t^2+5*t+13)*(t^4+7*t^3+20*t^2+19*t+1)^3/t;
    assert #Roots(Numerator(jt-0)) eq 0 and #Roots(Numerator(jt-1728)) eq 0;
    E:=EllipticCurveFromjInvariant(jt);
    assert IsSquare(Discriminant(E)/(t*(t^2+6*t+13)));
    R<x>:=PolynomialRing(Rationals());
    E:=EllipticCurve(x*(x^2-6*x+13));
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 2;
end procedure;

procedure VerifyRemark63()
    R<x>:=PolynomialRing(Rationals());
    F:=NumberField(x^3+x^2-2*x-1);
    E:=EllipticCurve("147b1");
    assert Q3inftyTorsion(E) eq [2,26] and #TorsionSubgroup(ChangeRing(E,F)) eq 13;
end procedure;

procedure VerifyLemma64()
    assert #{jInvariant(E): E in EllipticCurves(CremonaDatabase(),162,2)} eq 4;
    for E in EllipticCurves(CremonaDatabase(),162,2) do
        assert Q3inftyTorsion(E) eq [6,42];
    end for;
end procedure;

procedure VerifyLemma65()
    assert Q3inftyTorsion("49a1") eq [2,14];
    assert Q3inftyTorsion("49a2") eq [2,14];
end procedure;

procedure VerifyLemma68()
    assert #{jInvariant(E): E in EllipticCurves(CremonaDatabase(),50,1)} eq 4;
    assert Q3inftyTorsion("50a1") eq [6,6];
    assert Q3inftyTorsion("50a2") eq [6,6];
    assert Q3inftyTorsion("50a3") eq [6,30];
    assert Q3inftyTorsion("50a4") eq [6,30];
end procedure;

procedure VerifyProposition69()
    F<t>:=FunctionField(Rationals());
    jt:=(t^2+10*t+5)^3/t;
    assert #Roots(Numerator(jt-0)) eq 0 and #Roots(Numerator(jt-1728)) eq 0;
    E:=EllipticCurveFromjInvariant(jt);
    assert IsSquare(Discriminant(E)/(t*(t^2+22*t+125)));
    R<x>:=PolynomialRing(Rationals());
    E:=EllipticCurve(x*(x^2-22*x+125));
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 2;
end procedure;    
    
procedure VerifyLemma611()
    G4:=GL(2,Integers(4));
    B:=GL2Borel(Integers(4));
    B1:=GL2Borel1(Integers(4));
    G2,pi:=ChangeRing(G4,Integers(2));
    S:=[G`subgroup:G in Subgroups(G4)|#pi(G`subgroup) le 2];
    for G in S do
        N:=GenS3Divisor(G);
        if IsConjugateToSubgroup(G4,N,B1) then
            assert #pi(G) eq 1 or IsConjugateToSubgroup(G4,G,B);
        end if;
    end for;
end procedure;

procedure VerifyLemma612()
    assert Q3inftyTorsion("27a3") eq [18,18];
    F<t>:=FunctionField(Rationals());
    jt:=t^3*(t^3-24)^3/(t^3-27);
    assert #Roots(Numerator(jt-1728)) eq 0;
    R<x>:=PolynomialRing(Rationals());
    E:=EllipticCurve((x+3)*(x^2-3*x+9));
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 2;
end procedure;

procedure VerifyLemma613()
    print "This will take a minute...";
    F<t>:=FunctionField(Rationals());
    G9:=GL(2,Integers(9));
    S9:=SL(2,Integers(9));
    G3,pi:=ChangeRing(G9,Integers(3));
    B3:=GL2Borel(3);
    C3:=GL2SplitCartan(3);
    B9:=GL2Borel(9);
    T9:=GL2Borel1(9);
    H1:=sub<G9|[1,0,0,2],[2,0,0,1],[1,3,0,1],[1,1,6,1]>;
    H2:=sub<G9|[1,0,0,2],[2,0,0,1],[1,3,0,1],[1,1,3,1]>;
    H3:=sub<H2|[1,0,0,2],[2,0,0,1],[1,3,0,1]>;
    assert H3 eq H1 meet H2;
    assert IsConjugate(S9,H1 meet S9, CPGroups["9C0"]) and IsConjugate(S9,H2 meet S9, CPGroups["9A1"]) and IsConjugate(S9,H3 meet S9, CPGroups["9C1"]);
    S:=[G`subgroup:G in Subgroups(G9) | FullDeterminantMap(G`subgroup) and ContainsComplexConjugation(G`subgroup)];
    S:=[G:G in S | IsConjugateToSubgroup(G3,pi(G),B3) and IsConjugateToSubgroup(G9,GenS3Divisor(G),T9) and not IsConjugateToSubgroup(G3,pi(G),C3) and not IsConjugateToSubgroup(G9,G,B9) ];
    for G in S do assert IsConjugateToSubgroup(G9,G,H1) or IsConjugateToSubgroup(G9,G,H2); end for;
    assert RationalCuspCount(H2) eq 2 and RationalCuspCount(H3) eq 2;
    // Verify that Gamma0(3,9) is isomorphic to 27a1 as claimed by constructing it as a fiber product of maps from Gamma0(3,3) and Gamma0(9) over Gamma0(3)
    assert IsConjugate(G3,C3,SZGroups[[3,7]][2]) and IsConjugate(G9,B9,SZGroups[[9,292]][2]) and IsConjugate(G3,B3,SZGroups[[3,9]][2]);
    h1:=SZGroups[[3,7]][3];  h2:=SZGroups[[9,292]][3];  g:=SZGroups[[3,9]][3];
    f1:=-3*(t+1)^3/(t^2-t+1);  f2:=t*(t^2+9*t+27);
    assert h1 eq Evaluate(g,f1) and h2 eq Evaluate(g,f2);
    R<x,y>:=PolynomialRing(Rationals(),2);
    h:=Numerator(Evaluate(f1,x)-Evaluate(f2,y));
    assert IsIrreducible(h);
    C:=ProjectiveClosure(Curve(AffineSpace(R),h));
    assert Genus(C) eq 1;
    P:=[P:P in RationalPoints(C:Bound:=100)| not IsSingular(P)][1];
    E:=EllipticCurve(C,P);
    assert CremonaReference(E) eq "27a1";
    // Now show that 27a3 has only 3 rational points, 2 of which are cuspidal
    assert Rank(EllipticCurve("27a3")) eq 0 and #TorsionSubgroup(EllipticCurve("27a3")) eq 3;
    // this leaves H1 as the only possibility
    assert IsConjugate(G9,H1,SZGroups[[9,299]][2]);
    assert (t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3 eq SZGroups[[9,299]][3];
end procedure;

procedure VerifyExample614()
    F<t>:=FunctionField(Rationals());
    jt:=(t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3;
    E:=EllipticCurve("722a1");
    assert Evaluate(jt,-2) eq jInvariant(E);
    R<x>:=PolynomialRing(Rationals());
    K1:=NumberField(x^3 - x^2 - 6*x + 7);
    assert Discriminant(RingOfIntegers(K1)) eq 361;
    K2:=NumberField(x^3 - x^2 - 6*x - 12);
    assert Discriminant(RingOfIntegers(K2)) eq -1083;
    F:=Compositum(K1,K2);
    assert IsDivisibleBy(Invariants(TorsionSubgroup(ChangeRing(E,F)))[1], 9);
end procedure;

procedure VerifyProposition616()
    print "This will take a few minutes...";
    F<t,s>:=PolynomialRing(Rationals(),2);
    // Case 1 where E(Q)[2] is trivial
    jt := -4*(t^2-3)^3*(t^2-8*t-11)/(t+1)^4; // rational map corresponding to Ht
    js:=(s+27)*(s+3)^3/s;
    f:=Numerator(jt-js);
    assert IsIrreducible(f);
    C1:=ProjectiveClosure(Curve(AffineSpace(F),f));
    assert Genus(C1) eq 1;
    E,pi:=EllipticCurve(C1,C1![-1,0]);
    assert CremonaReference(E) eq "48a3";
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 8;
    // check that X_G has 4 rational cusps
    G12:=GL(2,Integers(12));
    G3,red3:=ChangeRing(G12,Integers(3));
    G4,red4:=ChangeRing(G12,Integers(4));
    Hs:=GL2Borel(3);
    Ht:=sub<G4|[3,0,0,3],[3,1,0,1],[0,3,1,3]>;
    G:=sub<G12|Kernel(red3),Inverse(red3)(Hs)> meet sub<G12|Kernel(red4),Inverse(red4)(Ht)>;
    assert RationalCuspCount(G) eq 4;
    // show that we have 4 points on C1 that map to non-cuspidal rational points on E
    S:={C1![-5,-36],C1![7,-81/4],C1![-5/4,-81/4],C1![-1/2,-36]};
    assert #{pi(P):P in S} eq 4;
    for P in S do assert Evaluate(js,[P[1],P[2]]) eq Evaluate(jt,[P[1],P[2]]); end for;
    // show that these 4 non-cuspidal points correspond to just 2 Qbar-isomorphism clases
    assert {Evaluate(js,[P[1],P[2]]) : P in S} eq {-35937/4,109503/64};
    // pick representatives for the two classes
    E1:=EllipticCurve("162a1");  E2:=EllipticCurve("162a2");
    assert {jInvariant(E1),jInvariant(E2)} eq {-35937/4,109503/64};
    // Check that these elliptic curves do not have any points of order 36 over Q(3^infty)
    assert not IsDivisibleBy(Q3inftyTorsion(E1)[2],36) and not IsDivisibleBy(Q3inftyTorsion(E1)[2],36);

    // Case 2 where E(Q)[2] is non-trivial
    js:=(s+256)^3/s^2;
    jt:=(t+3)*(t^2-3*t+9)*(t^3+3)^3 / t^3;
    f:=Numerator(jt-js);
    assert IsIrreducible(f);
    C2:=ProjectiveClosure(Curve(AffineSpace(F), f));
    assert Genus(C2) eq 2;
    b,C3,pi1 := IsHyperelliptic(C2);
    assert b;
    R<x>:=PolynomialRing(Rationals());
    C4:=HyperellipticCurve(x^6-34*x^3+1);
    b,pi2:=IsIsomorphic(C3,C4);
    assert b;
    assert RankBound(Jacobian(C4)) eq 0;
    S:=Chabauty0(Jacobian(C4));
    assert #S eq 6;
    // Verify that 4 of these 6 points are rational cusps on our smooth projective model C4 for X_H (fiber product of  the maps js and jt over X(1))
    // Here H is the intersection in GL(2,Z/18Z) of the full inverse images of our level 2 and level 9 subgroups
    G18:=GL(2,Integers(18));
    G2,red2:=ChangeRing(G18,Integers(2));
    G9,red9:=ChangeRing(G18,Integers(9));
    Hs:=sub<G2|[1,1,0,1]>;
    Ht:=sub<G9|[1,0,0,2],[2,0,0,1],[1,3,0,1],[1,1,6,1]>;
    G:=sub<G18|Kernel(red2),Inverse(red2)(Hs)> meet sub<G18|Kernel(red9),Inverse(red9)(Ht)>;
    assert RationalCuspCount(G) eq 4;
    // Now verity that we can account for the 2 non-cuspidal points on our smooth model for X_H
    P1:=C2![3,-16];  P2:=C2![-3,-256];
    // verify that P1 and P2 map to distinct rational points 
    assert pi2(pi1(P1)) ne pi2(pi1(P2)) and pi2(pi1(P1)) in S and pi2(pi1(P2)) in S;
    // verify that P1 and P2 correspond to non-cuspidal points on X_H by showing that the give points on both our genus 0 curves
    j1:=Evaluate(jt,[P1[1],P1[2]]);  j2:=Evaluate(jt,[P2[1],P2[2]]);
    assert j1 eq Evaluate(js,[P1[1],P1[2]]) and j2 eq Evaluate(js,[P1[1],P2[2]]);
    assert {j1,j2} eq {0,54000};
    // check that the corresponding elliptic curves do not have any points of order 36 over Q(3^infty)
    E1:=EllipticCurveFromjInvariant(j1);  E2:=EllipticCurveFromjInvariant(j2);
    assert not IsDivisibleBy(Q3inftyTorsion(E1)[2],36) and not IsDivisibleBy(Q3inftyTorsion(E1)[2],36);
end procedure;

// List of invariants of possible torsion structures of E(Q(3^infty)) for E/Q listed in Theorem 1.8
Tlist:=[[2,2],[2,4],[2,8],[2,10],[2,14],[2,16],[2,26],[4,4],[4,8],[4,16],[4,28],[6,6],[6,12],[6,18],[6,30],[6,42],[8,8],[12,12],[14,14],[18,18]];

// List of 33 maximal groups cited in the Proof of Theorem 7.1 -- these are verifed by VerifyTheorem71Groups() below
Glist:=AssociativeArray();
Glist[[2,2]]:=[sub<GL(2,Integers(2))|[1,1,0,1],[0,1,1,0]>];
Glist[[2,4]]:=[sub<GL(2,Integers(4))|[1,0,0,3],[3,0,0,1],[1,1,0,1]>];
Glist[[2,8]]:=[sub<GL(2,Integers(8))|[1,1,0,1],[1,0,0,5],[1,0,0,3],[3,0,0,1],[5,0,0,1]>];
Glist[[2,10]]:=[sub<GL(2,Integers(5))|[1,1,0,1],[1,0,0,2],[4,0,0,1]>];
Glist[[2,14]]:=[sub<GL(2,Integers(7))|[1,1,0,1],[1,0,0,3],[3,0,0,1]>];
Glist[[2,16]]:=[
    sub<GL(2,Integers(16))|[1,1,0,1],[9,0,0,1],[7,0,0,1],[1,0,0,7],[1,0,0,3]>,
    sub<GL(2,Integers(16))|[9,0,0,1],[1,2,0,1],[7,0,0,1],[1,0,0,7],[1,0,0,3],[3,1,8,1]>];
Glist[[2,26]]:=[sub<GL(2,Integers(13))|[1,1,0,1],[1,0,0,2],[4,0,0,1]>];
Glist[[4,4]]:=[
    sub<GL(2,Integers(2))|>,
    sub<GL(2,Integers(4))|[1,1,0,3],[3,0,0,3],[1,2,0,1]>,
    sub<GL(2,Integers(4))|[1,1,0,3],[3,0,0,3],[1,0,1,3]>];
Glist[[4,8]]:=[
    sub<GL(2,Integers(4))|[1,0,0,3],[3,0,0,1],[1,2,0,1]>,
    sub<GL(2,Integers(8))|[1,1,0,7],[1,0,0,5],[1,2,0,1],[5,0,0,1],[3,0,0,3]>];
Glist[[4,16]]:=[
    sub<GL(2,Integers(8))|[1,0,0,5],[1,2,0,1],[7,0,0,1],[1,0,0,3]>,
    sub<GL(2,Integers(16))|[1,0,0,5],[9,0,0,1],[1,2,0,1],[7,0,0,7],[1,1,0,15]>];
Glist[[4,28]]:=[
    sub<GL(2,Integers(14))|[1,2,0,1],[1,0,0,3],[3,0,0,1]>,
    sub<GL(2,Integers(28))|[1,0,0,5],[13,0,0,1],[3,0,0,3],[1,1,0,15],[1,0,7,3]>,
    sub<GL(2,Integers(28))|[1,7,0,15],[1,0,0,5],[1,2,0,1],[13,0,0,1],[3,0,0,3]>];
Glist[[6,6]]:=[sub<GL(2,Integers(3))|[1,1,0,1],[1,0,0,2],[2,0,0,1]>];
Glist[[6,12]]:=[sub<GL(2,Integers(12))|[1,1,0,1],[1,0,0,5],[7,0,0,1],[1,0,0,7],[5,0,0,1]>];
Glist[[6,18]]:=[
    sub<GL(2,Integers(9))|[1,1,0,1],[1,0,0,2],[2,0,0,1]>,
    sub<GL(2,Integers(9))|[1,0,0,2],[2,0,0,1],[1,3,0,1],[1,1,3,2]>,
    sub<GL(2,Integers(9))|[1,0,0,2],[2,0,0,1],[1,3,0,1],[1,1,3,1]>];
Glist[[6,30]]:=[sub<GL(2,Integers(15))|[1,1,0,1],[1,0,0,2],[11,0,0,1],[4,0,0,1],[1,0,0,11]>];
Glist[[6,42]]:=[sub<GL(2,Integers(21))|[1,1,0,1],[1,0,0,13],[1,0,0,2],[2,0,0,1],[13,0,0,1]>];
Glist[[8,8]]:=[sub<GL(2,Integers(4))|[1,0,0,3],[3,0,0,1]>];
Glist[[12,12]]:=[
    sub<GL(2,Integers(6))|[1,0,0,5],[1,2,0,1],[5,0,0,1]>,
    sub<GL(2,Integers(12))|[1,0,0,5],[1,2,0,1],[7,0,0,7],[1,3,0,7],[5,0,0,1]>,
    sub<GL(2,Integers(12))|[1,1,0,7],[1,0,0,5],[7,0,0,7],[1,0,3,7],[5,0,0,1]>];
Glist[[14,14]]:=[
    sub<GL(2,Integers(7))|[1,0,0,3],[3,0,0,1]>,
    sub<GL(2,Integers(7))|[1,0,0,2],[0,1,1,0],[3,0,0,3]>];
Glist[[18,18]]:=[
    sub<GL(2,Integers(3))|[1,0,0,2],[2,0,0,1]>,
    sub<GL(2,Integers(9))|[1,3,0,1],[1,0,0,8],[1,1,3,8],[2,0,0,2]>];
    
// List of the rational functions that appear in Table 1 (some of which are constant functions)
F<t>:=FunctionField(Rationals());
Tparam:=AssociativeArray();
Tparam[[2,2]]:=<t>;
Tparam[[2,4]]:=<(t^2+16*t+16)^3/(t*(t+16))>;
Tparam[[2,8]]:=<(t^4-16*t^2+16)^3/(t^2*(t^2-16))>;
Tparam[[2,10]]:=<(t^4-12*t^3+14*t^2+12*t+1)^3/(t^5*(t^2-11*t-1))>;
Tparam[[2,14]]:=<(t^2+13*t+49)*(t^2+5*t+1)^3/t>;
Tparam[[2,16]]:=<(t^16-8*t^14+12*t^12+8*t^10-10*t^8+8*t^6+12*t^4-8*t^2+1)^3/(t^16*(t^4-6*t^2+1)*(t^2+1)^2*(t^2-1)^4)>;
Tparam[[2,26]]:=<(t^4-t^3+5*t^2+t+1)*(t^8-5*t^7+7*t^6-5*t^5+5*t^3+7*t^2+5*t+1)^3/(t^13*(t^2-3*t-1))>;
Tparam[[4,4]]:=<(t^2+192)^3/(t^2-64)^2, -16*(t^4-14*t^2+1)^3/(t^2*(t^2+1)^4), -4*(t^2+2*t-2)^3*(t^2+10*t-2)/t^4>;
Tparam[[4,8]]:=<16*(t^4+4*t^3+20*t^2+32*t+16)^3/(t^4*(t+1)^2*(t+2)^4), -4*(t^8-60*t^6+134*t^4-60*t^2+1)^3/(t^2*(t^2-1)^2*(t^2+1)^8)>;
Tparam[[4,16]]:=<(t^16-8*t^14+12*t^12+8*t^10+230*t^8+8*t^6+12*t^4-8*t^2+1)^3/(t^8*(t^2-1)^8*(t^2+1)^4*(t^4-6*t^2+1)^2)>;
Tparam[[4,28]]:=<F!351/4, F!-38575685889/16384>;
Tparam[[6,6]]:=<(t+27)*(t+3)^3/t>;
Tparam[[6,12]]:=<(t^2-3)^3*(t^6-9*t^4+3*t^2-3)^3/(t^4*(t^2-9)*(t^2-1)^3)>;
Tparam[[6,18]]:=<(t+3)^3*(t^3+9*t^2+27*t+3)^3/(t*(t^2+9*t+27)), (t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3>;
Tparam[[6,30]]:=<F!-121945/32, F!46969655/32768>;
Tparam[[6,42]]:=<F!3375/2, F!-140625/8, F!-1159088625/2097152, F!-189613868625/128>;
Tparam[[8,8]]:=<(t^8+224*t^4+256)^3/(t^4*(t^4-16)^4)>;
Tparam[[12,12]]:=<(t^2+3)^3*(t^6-15*t^4+75*t^2+3)^3/(t^2*(t^2-9)^2*(t^2-1)^6),F!-35937/4, F!109503/64>;
Tparam[[14,14]]:=<F!2268945/128>;
Tparam[[18,18]]:=<27*t^3*(8-t^3)^3/(t^3+1)^3, 432*t*(t^2-9)*(t^2+3)^3*(t^3-9*t+12)^3*(t^3+9*t^2+27*t+3)^3*(5*t^3-9*t^2-9*t-3)^3/((t^3-3*t^2-9*t+3)^9*(t^3+3*t^2-9*t-3)^3)>;


procedure VerifyTheorem71Groups()
    print "Computing maximal subgroups G for each T -- will take several hours (about six hours on a 2.6GHz Intel Xeon E5-2670)...";
    Hlist:=AssociativeArray();
    for T in Tlist do
        Hlist[T] := MaximalImagesWithGenS3Torsion(T);
        assert #Hlist[T] eq #Glist[T];
        for H in Hlist[T] do assert #[G:G in Glist[T]|BaseRing(G) eq BaseRing(H) and IsConjugate(GL(2,BaseRing(G)),G,H)] eq 1; end for;
    end for;
    assert Keys(Hlist) eq Keys(Glist);
end procedure;
    
procedure VerifyTheorem71Maps()
    print "Verifying maps to X(1), this will take a few minutes...";
    assert #Glist[[2,4]] eq 1 and IsConjugate(GL(2,Integers(4)),Glist[[2,4]][1],SZGroups[[4,46]][2]) and Tparam[[2,4]][1] eq SZGroups[[4,46]][3];
    assert #Glist[[2,8]] eq 1 and IsConjugate(GL(2,Integers(8)),Glist[[2,8]][1],SZGroups[[8,2142]][2]) and Tparam[[2,8]][1] eq SZGroups[[8,2142]][3];
    assert #Glist[[2,10]] eq 1 and IsConjugate(GL(2,Integers(5)),Glist[[2,10]][1],SZGroups[[5,38]][2]) and Tparam[[2,10]][1] eq SZGroups[[5,38]][3];
    assert #Glist[[2,14]] eq 1 and IsConjugate(GL(2,Integers(7)),Glist[[2,14]][1],SZGroups[[7,76]][2]) and Tparam[[2,14]][1] eq SZGroups[[7,76]][3];
    assert #Glist[[2,16]] eq 2
        and IsConjugate(GL(2,Integers(16)),Glist[[2,16]][1],SZGroups[[16,25902]][2]) and Tparam[[2,16]][1] eq SZGroups[[16,25902]][3]
        and IsConjugate(GL(2,Integers(16)),Glist[[2,16]][2],RZBGroups[["H335"]]);
    assert #Glist[[2,26]] eq 1 and IsConjugate(GL(2,Integers(13)),Glist[[2,26]][1],SZGroups[[13,204]][2]) and Tparam[[2,26]][1] eq SZGroups[[13,204]][3];
    assert #Glist[[4,4]] eq 3
        and IsConjugate(GL(2,Integers(2)),Glist[[4,4]][1],SZGroups[[2,0]][2]) and Tparam[[4,4]][1] eq SZGroups[[2,0]][3]
        and IsConjugate(GL(2,Integers(4)),Glist[[4,4]][2],SZGroups[[4,32]][2]) and Tparam[[4,4]][2] eq SZGroups[[4,32]][3]
        and IsConjugate(GL(2,Integers(4)),Glist[[4,4]][3],SZGroups[[4,42]][2]) and Tparam[[4,4]][3] eq SZGroups[[4,42]][3];
    assert #Glist[[4,8]] eq 2
        and IsConjugate(GL(2,Integers(4)),Glist[[4,8]][1],SZGroups[[4,30]][2]) and Tparam[[4,8]][1] eq SZGroups[[4,30]][3]
        and IsConjugate(GL(2,Integers(8)),Glist[[4,8]][2],SZGroups[[8,1942]][2]) and Tparam[[4,8]][2] eq SZGroups[[8,1942]][3];
    assert #Glist[[4,16]] eq 2
        and IsConjugate(GL(2,Integers(8)),Glist[[4,16]][1],SZGroups[[8,1545]][2]) and Tparam[[4,16]][1] eq SZGroups[[8,1545]][3]
        and IsConjugate(GL(2,Integers(16)),Glist[[4,16]][2],RZBGroups[["H478"]]);

    assert #Glist[[4,28]] eq 3
        and IsConjugate(GL(2,Integers(2)),ChangeRing(Glist[[4,28]][1],Integers(2)),GL2SplitCartan(2))
        and IsConjugate(GL(2,Integers(7)),ChangeRing(Glist[[4,28]][1],Integers(7)),GL2Borel(7))
        and IsConjugateToSubgroup(GL(2,Integers(28)),Glist[[4,28]][3],GL2Borel(28));
    // for Glist[[4,28]][1] and Glist[[4,28]][3], if E has G as its Galois image then either E admits 2 2-isogenies and a 7-isogeny and it is 2-isogenous
    // to a curve E' that admits a 28-isogeny, or E admits a 28-isogeny, and neither can happen (no E/Q admits a rational 28-isogeny).
    // for Glist[[4,28]][2] we use map to X_H where H=[4,55]x[7,76]
    G28:=GL(2,Integers(28));
    G4,pi4:=ChangeRing(G28,Integers(4));  G7,pi7:=ChangeRing(G28,Integers(7));
    H:=ConjugateIntersection(G28,sub<G28|Kernel(pi4),Inverse(pi4)(SZGroups[[4,55]][2])>,sub<G28|Kernel(pi7),Inverse(pi7)(SZGroups[[7,76]][2])>);
    assert IsConjugate(SL(2,Integers(28)),H meet SL(2,Integers(28)),CPGroups["28A2"]);
    assert IsConjugateToSubgroup(G28,Glist[[4,28]][2],H);
    f:=SZGroups[[4,55]][3];  g:=SZGroups[[7,76]][3];
    R<x,y>:=PolynomialRing(Rationals(),2);
    h:=Numerator(Evaluate(f,x)-Evaluate(g,y));
    assert IsIrreducible(h);
    C:=ProjectiveClosure(Curve(AffineSpace(R),h));
    assert Genus(C) eq 2;
    b,XH:=IsHyperelliptic(C);
    R<z>:=PolynomialRing(Rationals());
    X:=HyperellipticCurve(z^6-2*z^5-4*z^4-4*z^3-4*z^2-2*z+1);
    assert b and IsIsomorphic(X,XH);
    J:=Jacobian(X);
    assert RankBounds(J) eq 1;
    P:=[P : P in  RationalPoints(J:Bound:=100) | Order(P) eq 0][1];
    assert #Chabauty(P) eq 4;   // 4 rational points on X_H
    assert RationalCuspCount(H) eq 2; // 2 are cusps
    // for all EQ with E(Q(3^infty))_tor eq [4,28] we know that j(E) lies in the image of X_G --> X_H --> X(1), and comes from a non-cuspidal rational point on X_H
    // there can thus be at most 2, and we know verify that there are in fact 2
    assert #[j:j in {351/4,-38575685889/16384} | Q3inftyTorsion(EllipticCurveFromjInvariant(j)) eq [4,28]] eq 2;
    
    assert #Glist[[6,6]] eq 1 and IsConjugate(GL(2,Integers(3)),Glist[[6,6]][1],SZGroups[[3,9]][2]) and Tparam[[6,6]][1] eq SZGroups[[3,9]][3];
    assert #Glist[[6,12]] eq 1 and IsConjugate(GL(2,Integers(12)),Glist[[6,12]][1],GL2Borel(12));   // j-map for X0(12) taken from Table 3 of Lozano
    assert #Glist[[6,18]] eq 3
        and IsConjugate(GL(2,Integers(9)),Glist[[6,18]][1],SZGroups[[9,292]][2]) and Tparam[[6,18]][1] eq SZGroups[[9,292]][3]
        and IsConjugate(GL(2,Integers(9)),Glist[[6,18]][2],SZGroups[[9,299]][2]) and Tparam[[6,18]][2] eq SZGroups[[9,299]][3]
        and IsConjugate(GL(2,Integers(9)),Glist[[6,18]][3],sub<GL(2,Integers(9))|[1,1,3,1],[2,0,0,1],[2,0,0,2]>); // H3 from Lemma 5.21, verified above
    assert #Glist[[6,30]] eq 1 and IsConjugateToSubgroup(GL(2,Integers(15)),Glist[[6,30]][1],GL2Borel(15));
    // The group G is the direct product [3,9]x[5,38],which is conjugate to an index-2  subgroup of the Borel in GL(2,Z/15Z)
    // we don't need to compute a model for X_G, we can just check the 4 j-invariants that admit a rational 15-isogeny (represented by isogeny classes 50a)
    assert {jInvariant(E):E in EllipticCurves(CremonaDatabase(),50,1)|Q3inftyTorsion(E) eq [6,30]} eq {-121945/32, 46969655/32768};
    assert #Glist[[6,42]] eq 1 and IsConjugate(GL(2,Integers(21)),Glist[[6,42]][1],GL2Borel(21));
    // Every j(E) corresponding to a rational point on X_0(21) must have E(Q(3^infty))_tor = [6,42], and we know there are 4, but we may well verify
    assert #[j:j in {3375/2,-140625/8,-1159088625/2097152,-189613868625/128} | Q3inftyTorsion(EllipticCurveFromjInvariant(j)) eq [6,42]] eq 4;
    assert #Glist[[8,8]] eq 1 and IsConjugate(GL(2,Integers(4)),Glist[[8,8]][1],SZGroups[[4,15]][2]) and Tparam[[8,8]][1] eq SZGroups[[4,15]][3];
    
    assert #Glist[[12,12]] eq 3;
    // the first group has level 6 and genus 0 -- H=[2,0]x[3,9] and we construct the map to X(1) as a fiber product
    H:=Glist[[12,12]][1];
    G6:=GL(2,Integers(6));  G2,pi2:=ChangeRing(G6,Integers(2));  G3,pi3:=ChangeRing(G6,Integers(3));
    assert IsConjugate(G6,H,ConjugateIntersection(G6,sub<G6|Kernel(pi2),Inverse(pi2)(SZGroups[[2,0]][2])>,sub<G6|Kernel(pi3),Inverse(pi3)(SZGroups[[3,9]][2])>));
    f:=SZGroups[[2,0]][3];  g:=SZGroups[[3,9]][3];
    R<x,y>:=PolynomialRing(Rationals(),2);
    F<t>:=FunctionField(Rationals());
    h:=Numerator(Evaluate(f,x)-Evaluate(g,y));
    assert IsIrreducible(h);
    C:=ProjectiveClosure(Curve(AffineSpace(R),h));
    assert Genus(C) eq 0;
    Co,phi1:=Conic(C);
    assert HasPoint(Co);
    P1:=Curve(ProjectiveSpace(Rationals(),1));
    phi2:=Parametrization(Co,P1);
    phi3:=Expand(phi2*Inverse(phi1));
    phi3polys:=DefiningPolynomials(phi3);
    h:=Evaluate(f,phi3polys[1]/phi3polys[3]);
    j:=Evaluate(h,[t,1]);
    // Due to randomness in Magma, the map we get will vary, so just do a quick sanity check to convince ourselves it is equivalent to Tparam[[12,12]][1]
    f:=Tparam[[12,12]][1];
    assert Degree(f) eq Degree(j);
    for i:=1 to 100 do p:=Random(1000000); if Evaluate(Denominator(f),p) ne 0 then assert #Roots(Numerator(j-Evaluate(f,p))) gt 0; end if; end for;
    for i:=1 to 100 do p:=Random(1000000); if Evaluate(Denominator(j),p) ne 0 then assert #Roots(Numerator(f-Evaluate(j,p))) gt 0; end if; end for;
    // the second group has level 12  and genus 1 -- H=[4,32]x[3,39] and X_H is isomorphic to 48a1 with 4 rational points, all cuspidal
    H:=Glist[[12,12]][2];
    G12:=GL(2,Integers(12));  G4,pi4:=ChangeRing(G12,Integers(4));  G3,pi3:=ChangeRing(G12,Integers(3));
    assert IsConjugate(G12,H,ConjugateIntersection(G12,sub<G12|Kernel(pi4),Inverse(pi4)(SZGroups[[4,32]][2])>,sub<G12|Kernel(pi3),Inverse(pi3)(SZGroups[[3,9]][2])>));
    f:=SZGroups[[4,32]][3];  g:=SZGroups[[3,9]][3];
    R<x,y>:=PolynomialRing(Rationals(),2);
    F<t>:=FunctionField(Rationals());
    h:=Numerator(Evaluate(f,x)-Evaluate(g,y));
    assert IsIrreducible(h);
    C:=ProjectiveClosure(Curve(AffineSpace(R),h));
    P:=[P:P in RationalPoints(C:Bound:=10)][1];
    E,phi:=EllipticCurve(C,P);
    assert CremonaReference(E) eq "48a1";
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 4;
    assert RationalCuspCount(H) eq 4;
    // the third group has level 12  and genus 1 -- H=[4,42]x[3,39] and X_H is isomorphic to 48a3 with 8 rational points, 4 cuspidal
    H:=Glist[[12,12]][3];
    G12:=GL(2,Integers(12));  G4,pi4:=ChangeRing(G12,Integers(4));  G3,pi3:=ChangeRing(G12,Integers(3));
    assert IsConjugate(G12,H,ConjugateIntersection(G12,sub<G12|Kernel(pi4),Inverse(pi4)(SZGroups[[4,42]][2])>,sub<G12|Kernel(pi3),Inverse(pi3)(SZGroups[[3,9]][2])>));
    f:=SZGroups[[4,42]][3];  g:=SZGroups[[3,9]][3];
    R<x,y>:=PolynomialRing(Rationals(),2);
    F<t>:=FunctionField(Rationals());
    h:=Numerator(Evaluate(f,x)-Evaluate(g,y));
    assert IsIrreducible(h);
    C:=ProjectiveClosure(Curve(AffineSpace(R),h));
    P:=[P:P in RationalPoints(C:Bound:=10)|not P in SingularPoints(C)][1];
    E,phi:=EllipticCurve(C,P);
    assert CremonaReference(E) eq "48a3";
    assert Rank(E) eq 0 and #TorsionSubgroup(E) eq 8;
    assert RationalCuspCount(H) eq 4;
    S:=[P:P in RationalPoints(C:Bound:=1000)| P[3] ne 0 and Evaluate(Denominator(f),P[1]/P[3]) ne 0 and Evaluate(Denominator(g),P[2]/P[3]) ne 0];   // rational affine points where we can evaluate both f and g give us j-invariants
    assert #{phi(P):P in S} eq 4;   // we get 4 distinct points on E, and none of tehm are cuspidal because they all map to E/Q with E(Q(3^infty))_tor = [12,12]
    assert {Evaluate(f,P[1]/P[3]):P in S} eq {-35937/4,109503/64};
    assert {Evaluate(g,P[2]/P[3]):P in S} eq {-35937/4,109503/64};
    assert #[j:j in {-35937/4,109503/64} | Q3inftyTorsion(EllipticCurveFromjInvariant(j)) eq [12,12]] eq 2;
    
    assert #Glist[[18,18]] eq 2
        and IsConjugate(GL(2,Integers(3)),Glist[[18,18]][1],SZGroups[[3,7]][2]) and Tparam[[18,18]][1] eq SZGroups[[3,7]][3]
        and IsConjugate(GL(2,Integers(9)),Glist[[18,18]][2],SZGroups[[9,268]][2]) and Tparam[[18,18]][2] eq SZGroups[[9,268]][3];
end procedure;

function Q3inftyTorsionFast(E)
    if (Type(E) eq RngIntElt or Type(E) eq FldRatElt) then
        if E eq 1728 then print "For j-invariant 1728 you must specify a particular curve E/Q"; assert false; end if;
        E:=EllipticCurveFromjInvariant(E);
    end if;
    if Type(E) ne CrvEll then E:=EllipticCurve(E); end if;
    assert BaseRing(E) eq Rationals();
    jE := jInvariant(E);
    if jE eq 1728 then return Q3inftyTorsion(E); end if;
    Tmax:=[1,1];
    for T in Keys(Tparam) do
        for f in Tparam[T] do
            if f-jE eq 0 or #Roots(Numerator(f-jE)) gt 0 then Tmax:=[LCM(Tmax[1],T[1]),LCM(Tmax[2],T[2])]; end if;
        end for;
    end for;
    assert Tmax in Keys(Tparam);
    return Tmax;
end function;

// Note that the two sanity checks below are slightly undermined by the fact that the code for Q3inftyTorsion actually relies on some of the parameterizations above to handle points of order 7 and 9.
// Unfortunately, it is not computationally feasible to run these tests without this special handling,   We  explicitly depend on Proposition 4.1, and Lemmas 5.13, 5.19, and 5.21, but nothing else.
// In particular, unless you set UseTheorem51 to true, the two checks below do not rely on any of our more delicate arguments (e.g. ruling out points of order 27 or 32)

// Verify that the parameterizations for torsion group T  produce a curve E/Q for which E(Q(3^infty))_tor contains T for n random values of T
// This is time consuming, so you will want to keep n small, say less than 10
function ParamCheck1(n:UseTheorem51:=false)
    printf "Sanity checking each parameterization at %o random j-invariants, this will take a while (possibly hours), be patient...\n", n;
    for T in Keys(Tparam) do
        print "Checking parameterization for T =", T;
        for f in Tparam[T] do
            k:=0; m:=0;
            for i:=1 to n do
                r:=Random(1,100)/Random(1,100);
                if Evaluate(Denominator(f),r) ne 0 then
                    E:=EllipticCurveFromjInvariant(Evaluate(f,r));
                    Etor:=Q3inftyTorsion(E:UseTheorem51:=UseTheorem51);
                    if not TorsionContains(Etor,T) then
                        printf "Verification for T = %o failed using function f = %o with t = %o yielding E = %o\n", T, f, r, E;
                        assert false;
                    end if;
                    k +:= 1;
                    if Etor eq T then m +:=1; end if;
                    if Degree(f) eq 0 then break; end if;
                end if;
            end for;
            printf "Torsion subgroup %o function %o passed %o tests successfully, with %o exact matches\n", T, f, k, m;
        end for;
    end for;
    return true;
end function;

// Given a conductor bound or list of conductors N, verify that for each elliptic curve  E of conductor <= maxN the j-invariant j(E) satisfies one of the parameterizations for the torsion group T = E(Q(3^infty))_tor
function ParamCheck2(N:UseTheorem51:=false)
    printf "Sanity checking parameterizations for curves of conductor up to %o\n", N;
    DB:=CremonaDatabase();
    if Type(N) ne SeqEnum then N := [11..N]; end if;
    for n in N do
        for E in EllipticCurves(DB,n) do
            T:=Q3inftyTorsion(E:UseTheorem51:=UseTheorem51);
            printf "%o %o\n", CremonaReference(E), T;
            assert Q3inftyTorsionFast(E) eq T;
            if jInvariant(E) eq 1728 then assert T in {[2,2],[4,4]}; end if;
        end for;
    end for;
    printf "Successfully checked E(Q(3^infty))_tor parameterizations for all elliptic curves E/Q with conductor in [%o...%o]\n", Min(N), Max(N);
    return true;
end function;

procedure VerifyAll()
    print "*** Verifying Example 1.9 ***";
    VerifyExample19();
    print "*** Example 1.9 Verified ***";
    print "*** Verifying Lemma 5.10 ***";
    VerifyLemma510();
    print "*** Lemma 5.10 Verified ***";
    print "*** Verifying Corollary 5.15 ***";
    VerifyCorollary515();
    print "*** Corollary 5.15 Verified ***";
    print "*** Verifying Lemma 5.16 ***";
    VerifyLemma516();
    print "*** Lemma 5.16 Verified ***";
    print "*** Verifying Lemma 5.17 ***";
    VerifyLemma517();
    print "*** Lemma 5.17 Verified ***";
    print "*** Verifying Lemma 5.18 ***";
    VerifyLemma518();
    print "*** Lemma 5.18 Verified ***";
    print "*** Verifying Lemma 5.19 ***";
    VerifyLemma519();
    print "*** Lemma 5.19 Verified ***";
    print "*** Verifying Example 5.20 ***";
    VerifyExample520();
    print "*** Example 5.20 Verified ***";
    print "*** Verifying Lemma 5.21 ***";
    VerifyLemma521();
    print "*** Lemma 5.21 Verified ***";
    print "*** Verifying Lemma 6.1 ***";
    VerifyLemma61();
    print "*** Lemma 6.1 Verified ***";
    print "*** Verifying Lemma 6.2 ***";
    VerifyProposition62();
    print "*** Proposition 6.2 Verified ***";
    print "*** Verifying Remark 6.3 ***";
    VerifyRemark63();
    print "*** Remark 6.3 Verified ***";
    print "*** Verifying Lemma 6.4 ***";
    VerifyLemma64();
    print "*** Lemma 6.4 Verified ***";
    print "*** Verifying Lemma 6.5 ***";
    VerifyLemma65();
    print "*** Lemma 6.5 Verified ***";
    print "*** Verifying Lemma 6.8 ***";
    VerifyLemma68();
    print "*** Lemma 6.8 Verified ***";
    print "*** Verifying Proposition 6.9 ***";
    VerifyProposition69();
    print "*** Proposition 6.9 Verified ***";
    print "*** Verifying Lemma 6.11 ***";
    VerifyLemma611();
    print "*** Lemma 6.11 Verified ***";
    print "*** Verifying Lemma 6.12 ***";
    VerifyLemma612();
    print "*** Lemma 6.12 Verified ***";
    print "*** Verifying Lemma 6.13 ***";
    VerifyLemma613();
    print "*** Lemma 6.13 Verified ***";
    print "*** Verifying Example 6.14 ***";
    VerifyExample614();
    print "*** Example 6.14 Verified ***";
    print "*** Verifying Proposition 6.16 ***";
    VerifyProposition616();
    print "*** Proposition 6.16 Verified ***";
    print "*** Verifying Theorem 7.1 ***";
    VerifyTheorem71Groups();
    VerifyTheorem71Maps();
    print "*** Theorem 7.1 Verified ***";
    print "*** Sanity checking parametrization in Theorem 7.1 ***";
    ParamCheck1(3);
    ParamCheck2(2500);
    print "*************************************************";
    print "*** Verification complete -- all tests passed ***";
    print "*************************************************";
end procedure;

print "Use VerifyAll() to verify all the Magma computations cited in the article -- this will take several hours.";
print "The function Q3inftyTorsion(E) computes E(Q(3^infty))_tor for any elliptic curve over E/Q, using the algorithm described in Section 5.5 of the paper.";
print "The function Q3inftyTorsionFast(E) is a faster version that uses the parametrizations of Theorem 7.1.";


