/*
    Magma code related to the paper "Torsion subgroups of rational elliptic curves over the compositum of all extensions of generalized $D_4$-type"
    by H. Daniels
    Copyright (c) 2017 the author.
    
    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)
	Much of this code has been taken from [5] and modified.


    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(D4^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).
*/

/*This code requires loading the database on the 727 arithmetically maximal subgroups containing −I that are in inside GL2(Z_2).
This data is available at

http://users.wfu.edu/rouseja/2adic/gl2data.gz 
*/ 

//load "path to gl2data";


// 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 4 -- avoids wasting time computing the splitting field K of f
// If the algorithm returns true then K is not of generalized D4-type, if it returns false then it may or may not be.
function NotGenD4Poly(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 4 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 D4-type via Theorem 3.2
function IsGenD4(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 NotGenD4Poly(DefiningPolynomial(G)) then return false; else return $$(GaloisGroup(G)); end if;
    end if;
    if Type(G) eq RngUPolElt then
        if NotGenD4Poly(G) then return false; else return $$(GaloisGroup(G)); end if;
    end if;
    return (4 mod Exponent(G)) eq 0 and IsNilpotent(G) and  NilpotencyClass(G) le 2;
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 D4-type
function GenD4Divisor(G)
    if IsGenD4(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 IsGenD4(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(D4^infty))
function GenD4Torsion(G)
    A := FixModule(GenD4Divisor(G));
    if #A eq 1 then A := [1,A[1]]; end if;
    return A;
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] ge 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(D4^infty))[n]
// is isomorphic to T.  The returned groups will each be reduced modulo m=GL2Level(G) (unless m=1)
function MaximalImagesWithGenD4Torsion(T)
    assert ValidTorsionSubgroup(T);
    assert #T eq 2;
    n := T[#T];
    A:=T;
    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|GenD4Torsion(H) eq A];
    t:=Cputime();
    printf "%o of %o candidate images G in GL(2,Z/%oZ) yielding torsion %o over Q(D4^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;
    

//The an associative array for the j-maps computed in Section 5.
//The keys are the associated torsion structures.
F<t> := FunctionField(Rationals());
pPrimMaps := AssociativeArray();
pPrimMaps[[ 1, 1 ]]:=<t>;
pPrimMaps[[ 1, 3 ]] := <(27*t^4 + 756*t^3 + 7290*t^2 + 26244*t + 19683)/t^3>;
pPrimMaps[[ 1, 5 ]] := <25 * (t^2+10*t+5)^3/t^5>;
pPrimMaps[[ 1, 7 ]] := <(t^2-t+1)^3*(t^6-11*t^5+30*t^4-15*t^3-10*t^2+5*t+1)^3*(t-1)^-7*t^-7*(t^3-8*t^2+5*t+1)^-1 >;
pPrimMaps[[ 1, 9 ]] := <(t^36 - 36*t^35 + 594*t^34 - 6000*t^33 + 41859*t^32 - 215892*t^31 + 860466*t^30 - 2733948*t^29 + 7083369*t^28 - 15215576*t^27 + 27417258*t^26 - 41737032*t^25 + 53766648*t^24 - 58258332*t^23 + 52042986*t^22 - 36346716*t^21 + 16629516*t^20 + 241380*t^19 - 9621660*t^18 + 11057688*t^17 - 7521093*t^16 + 2914848*t^15 + 178308*t^14 - 1208520*t^13 + 948843*t^12 - 375948*t^11 + 15714*t^10 + 78304*t^9 - 48717*t^8 + 11124*t^7 + 2262*t^6 - 2124*t^5 + 378*t^4 + 84*t^3 - 36*t^2 + 1)/(t^27 - 18*t^26 + 144*t^25 - 702*t^24 + 2376*t^23 - 5985*t^22 + 11688*t^21 - 18135*t^20 + 22662*t^19 - 22923*t^18 + 18711*t^17 - 12168*t^16 + 6132*t^15 - 2259*t^14 + 522*t^13 - 27*t^12 - 27*t^11 + 9*t^10 - t^9)>;
pPrimMaps[[ 1, 13]] := <(t^2 - t + 1)^3*( t^12 - 9*t^11 + 29*t^10 - 40*t^9 + 22*t^8 - 16*t^7 + 40*t^6 - 22*t^5 - 23*t^4 + 25*t^3 - 4*t^2 - 3*t + 1)^3 * (t - 1)^(-13)*t^(-13)*(t^3 - 4*t^2 + t +1)^-1>;
pPrimMaps[[ 3, 3 ]] := <(27*(t+1)^3*(t-3)^3)/t^3>;
pPrimMaps[[ 4, 4 ]] := <t^3/(t + 16)>;
pPrimMaps[[ 4, 8 ]] := <(t^6 + 48*t^4 + 768*t^2 + 4096)/t^2, (t^6 - 144*t^4 + 6912*t^2 - 110592)/(t^2 - 64)>;
pPrimMaps[[ 4, 16 ]] := <(t^12 - 48*t^10 + 816*t^8 - 5632*t^6 + 13056*t^4 - 12288*t^2 + 4096)/(t^4 - 16*t^2)>;
pPrimMaps[[ 4, 32 ]] := <(t^48 - 24*t^46 + 228*t^44 - 1064*t^42 + 2322*t^40 - 840*t^38 - 5676*t^36 + 8712*t^34 + 1455*t^32 - 16240*t^30 + 14664*t^28 + 9456*t^26 - 21892*t^24 + 9456*t^22 + 14664*t^20 - 16240*t^18 + 1455*t^16 + 8712*t^14 - 5676*t^12 - 840*t^10 + 2322*t^8 - 1064*t^6 + 228*t^4 - 24*t^2 + 1)/(t^32 - 8*t^30 + 12*t^28 + 8*t^26 - 26*t^24 + 8*t^22 + 12*t^20 - 8*t^18 + t^16)>;
pPrimMaps[[ 5, 5 ]] := <(t^30 + 30*t^29 + 465*t^28 + 4900*t^27 + 39150*t^26 + 252000*t^25 + 1360000*t^24 + 6337500*t^23 + 26090625*t^22 + 96593750*t^21 + 325846875*t^20 + 1010625000*t^19 + 2898984375*t^18 + 7721718750*t^17 + 19152890625*t^16 + 44323437500*t^15 + 95764453125*t^14 + 193042968750*t^13 + 362373046875*t^12 + 631640625000*t^11 + 1018271484375*t^10 + 1509277343750*t^9 + 2038330078125*t^8 + 2475585937500*t^7 + 2656250000000*t^6 + 2460937500000*t^5 + 1911621093750*t^4 + 1196289062500*t^3 + 567626953125*t^2 + 183105468750*t + 30517578125)/(t^25 + 25*t^24 + 325*t^23 + 2875*t^22 + 19250*t^21 + 103125*t^20 + 456875*t^19 + 1709375*t^18 + 5475000*t^17 + 15140625*t^16 + 36321875*t^15 + 75703125*t^14 + 136875000*t^13 + 213671875*t^12 + 285546875*t^11 + 322265625*t^10 + 300781250*t^9 + 224609375*t^8 + 126953125*t^7 + 48828125*t^6 + 9765625*t^5)>;
pPrimMaps[[ 8, 8 ]] := <(64*t^6 + 576*t^4 + 1728*t^2 + 1728)/(t^4 - 2*t^2 + 1), (t^6 - 48*t^4 + 768*t^2 - 4096)/t^2>;
pPrimMaps[[ 8, 16 ]] := <(256*t^12 - 768*t^10 + 1536*t^8 - 1792*t^6 + 1536*t^4 - 768*t^2 + 256)/(t^8 - 2*t^6 + t^4), (-4*t^24 + 720*t^22 - 44808*t^20 + 1057680*t^18 - 6090684*t^16 + 15714720*t^14 - 21291632*t^12 + 15714720*t^10 - 6090684*t^8 + 1057680*t^6 - 44808*t^4 + 720*t^2 - 4)/(t^22 + 6*t^20 + 13*t^18 + 8*t^16 - 14*t^14 - 28*t^12 - 14*t^10 + 8*t^8 + 13*t^6 + 6*t^4 + t^2)>;
pPrimMaps[[ 8, 32 ]] := <(256) * (t - 1/2)^-8 * t^-8 * (t + 1/2)^-8 * (t^2 - t - 1/4)^-2 * (t^2 + 1/4)^-4 * (t^2 + t - 1/4)^-2 * (t^16 - 2*t^14 + 3/4*t^12 + 1/8*t^10 + 115/128*t^8 + 1/128*t^6 + 3/1024*t^4 - 1/2048*t^2 + 1/65536)^3>;
pPrimMaps[[ 16, 16 ]] := <(16*t^24 + 672*t^20 + 9456*t^16 + 45248*t^12 + 9456*t^8 + 672*t^4 + 16)/(t^20 - 4*t^16 + 6*t^12 - 4*t^8 + t^4)>;






//List of j-invatiants of elliptic curves with CM
function HasCm(E)
	Cmjs := [-262537412640768000, -147197952000, -884736000, -12288000, -884736, -32768, -3375, 0, 1728,8000, 54000, 287496, 16581375];
    return jInvariant(E) in Cmjs;
end function;

//Representatives of each Q(D4^infty)-isomorphism class of CM Elliptic curve and their torsion structure over Q(D_^\infty) as listed in Table 3
CMTorsionArray := AssociativeArray();
	CMTorsionArray["26569a1"] := [ 1, 1 ];
	CMTorsionArray["4489a1"] := [ 1, 1 ];
	CMTorsionArray["1849a1"] := [ 1, 1 ];
	CMTorsionArray["27a4"] := [ 1, 3 ];
	CMTorsionArray["361a1"] := [ 1, 1 ];
	CMTorsionArray["121b1"] := [ 3, 3 ];
	CMTorsionArray["49a1"] := [ 8, 8 ];
	CMTorsionArray["27a1"] := [ 3, 3 ];
	CMTorsionArray["36a1"] := [ 8, 24 ];
	CMTorsionArray["108a2"] := [ 1, 3 ];
	CMTorsionArray["32a1"] := [ 8, 8 ];
	CMTorsionArray["256d1"] := [ 12, 24 ];
	CMTorsionArray["36a2"] := [ 8, 24 ];
	CMTorsionArray["32a4"] := [ 4, 8 ];
	CMTorsionArray["49a2"] := [ 8, 8 ];
    

//Given E/Q, determines if E(Q(D4^infy))[p^n]=E[p^n].
//If UseTheorem51 is true (default) use the j-maps from RZB when p^n in [8,16]
function HasD4inftyFullTorsion(E,p,n: UseRZBData:=true)
    assert IsPrime(p) and n ge 1;
    //Uses Lemma 3.10 and Proposition 4.1 to rule out p^n not in [2,3,4,5,8,16]
    if not p^n in [2,3,4,5,8,16] then return false; end if;
    //Use the j-Maps from RZB when p^n in [8,16]
    if UseRZBData and p^n eq 8 and jInvariant(E) ne 0 then
        F<t>:=FunctionField(Rationals());
        J8_1 := (t^6 - 48*t^4 + 768*t^2 - 4096)/t^2;
        J8_2 := (64*t^6 + 576*t^4 + 1728*t^2 + 1728)/(t^4 - 2*t^2 + 1);
        return #Roots(Numerator(J8_1-jInvariant(E))) + #Roots(Numerator(J8_2-jInvariant(E))) gt 0;
    end if;
    if UseRZBData and p^n eq 16 and jInvariant(E) ne 0 then
        F<t>:=FunctionField(Rationals());
        J16 := (16*t^24 + 672*t^20 + 9456*t^16 + 45248*t^12 + 9456*t^8 + 672*t^4 + 
        16)/(t^20 - 4*t^16 + 6*t^12 - 4*t^8 + t^4);
        return #Roots(Numerator(J16-jInvariant(E))) gt 0;
    end if;
    E := WeierstrassModel(E);
    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 IsGenD4(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 IsGenD4(L) then return false; end if;
    end for;
    return true;
end function;

//Given E/Q, determines if E(Q(D4^infy))(p) contains a point of order p^n.
//If UseTheorem51 is true (default) use the group Tmax in Theorem 5.1 to rule out and p^n that does not divide #Tmax
//If UsepPrimMaps is true (default) use the maps computed in Section 5
function HasD4inftyTorsionPoint(E,p,n:UseTmax:=true, UsepPrimMaps := true)
    assert IsPrime(p) and n ge 1;
    if not p in [2,3,5,7,13] then return false; end if;
    if p in [7,13] and n gt 1 then return false; end if;
    if UseTmax then
        if p^n gt 32 or p^n in [25,27] then return false; end if;
    end if;
    // Handle p^n eq 7, 9, 16, and 32 specially to keep compuation time low.
    if UsepPrimMaps and p^n eq 7 and jInvariant(E) ne 0 then
        F<t>:=FunctionField(Rationals());
        J7 := (t^2-t+1)^3*(t^6-11*t^5+30*t^4-15*t^3-10*t^2+5*t+1)^3*(t-1)^-7*t^-7*(t^3-8*t^2+5*t+1)^-1;
        return #Roots(Numerator(J7-jInvariant(E))) gt 0;
    end if;
    if UsepPrimMaps and p^n eq 9 and jInvariant(E) ne 0 then
        F<t>:=FunctionField(Rationals());
        J9 := pPrimMaps[[1,9]][1];
        return #Roots(Numerator(J9-jInvariant(E))) gt 0;
    end if;
    if UsepPrimMaps and p^n eq 16 and jInvariant(E) ne 0 then
        F<t>:=FunctionField(Rationals());
        J16_1 := pPrimMaps[[ 4, 16 ]][1];
        J16_2 := pPrimMaps[[ 8, 16 ]][1];
        J16_3 := pPrimMaps[[ 8, 16 ]][2];
        J16_4 := pPrimMaps[[ 16, 16 ]][1];
        return (#Roots(Numerator(J16_1-jInvariant(E))) + #Roots(Numerator(J16_2-jInvariant(E))) + #Roots(Numerator(J16_3-jInvariant(E))) + #Roots(Numerator(J16_4-jInvariant(E)))) gt 0;
    end if;
    if UsepPrimMaps and p^n eq 32 and jInvariant(E) ne 0 then
        F<t>:=FunctionField(Rationals());
        J32_1 := pPrimMaps[[ 4, 32 ]][1];
        J32_2:= pPrimMaps[[ 8, 32 ]][1];
        return (#Roots(Numerator(J32_1-jInvariant(E))) + #Roots(Numerator(J32_2-jInvariant(E)))) gt 0;
    end if;
    E := WeierstrassModel(E);
    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 IsGenD4(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 IsGenD4(L) then return true; end if;
    end for;
    return false;
end function;

//Given E/Q, determines  E(Q(D4^infy))(p)
//If UseTheorem51 is true (default) use the group Tmax in Theorem 5.1 to rule out and p^n that does not divide #Tmax
//If UsepPrimMaps is true (default) use the maps computed in Section 5
function pPrimaryD4inftyTorsion(E,p:UseTmax:=true,UseRZBData:=true,UsepPrimMaps:=true)
    E:=WeierstrassModel(IntegralModel(E));
    m:=0;
    while true do
        if not HasD4inftyFullTorsion(E,p,m+1: UseRZBData:=UseRZBData) then break; end if;
        m:=m+1;
    end while;
    n:=m;
    while true do
        if not HasD4inftyTorsionPoint(E,p,n+1:UseTmax:=UseTmax,UsepPrimMaps:=UsepPrimMaps) then break; end if;
        n:=n+1;
    end while;
    return [p^m,p^n];
end function;

//Given E/Q, determines E(Q(D4^infy))_tor
// If UseCmData is true, then it uses the precomputed CM torsion structures in Table 2.
function D4inftyTorsion(E:UseTmax:=true,UseRZBData:=true,UsepPrimMaps:=true, UseCMData:=true)
    // To compute these set all parameters to false.
    if UseCMData and HasCm(E) and jInvariant(E) ne 0 then
    	for k in Keys(CMTorsionArray) do
        	if IsQuadraticTwist(E, EllipticCurve(k)) then
            return CMTorsionArray[k];
            end if;
        end for;
    end if;
    if (Type(E) eq RngIntElt or Type(E) eq FldRatElt) then
        if E eq 0 then print "For j-invariant 0 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();
    for p in [2,3,5,7,13] do
        Tp:=pPrimaryD4inftyTorsion(E,p:UseTmax:=UseTmax,UsepPrimMaps:=UsepPrimMaps);
        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=2890
function D4inftyTorsionExamples(maxN)
    examples:=AssociativeArray();
    D:=CremonaDatabase();
    for N in [11..maxN] do
        for E in EllipticCurves(D,N) do
           T:=D4inftyTorsion(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;

//Given a subgroup G of GL(2,Z/mZ) lifts or reduces the group to level n
function ChangeLevel(G,n)
    I := BaseRing(G);
    if #I ge n then
    H := ChangeRing(G,Integers(n));
    end if;
    if not #I ge n then
    GL2n := GL(2,Integers(n));
    _,pi  := ChangeRing(GL(2,Integers(n)),I);
    H := sub<GL2n | Inverse(pi)(G),Kernel(pi) >;
    end if;
    return H;
end function;


//Given a subgroup G of GL(2,m) lifts or reduces the group to level n and returns [a,b] with a|b|n
// specifying the subgroup of Z/nZ x Z/nZ fixed by Gal(Q(E[n]) cap Q(D4^infty))
function TwoPowerLevelStructure(G,n)
    H := ChangeLevel(G,n);
    return GenD4Torsion(H);
end function;


function MatrixGroupTranspose(H)
	G := GL(2,BaseRing(H));
    return sub<G| [Transpose(A) : A in Generators(H) ]>;
end function;


// 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 VerifyTable1()
    print "This will take several minutes...";
    A:=D4inftyTorsionExamples(338);
    E:=EllipticCurve("2890d1");
    //E:=MinimalTwist(EllipticCurveFromjInvariant(-60698457/40960));
    A[D4inftyTorsion(E)]  := [CremonaReference(E)];
    assert #A eq 24;
    assert A[[1,1]][1] eq "26b2";
    assert A[[1,3]][1] eq "19a2";
    assert A[[1,5]][1] eq "11a2";
    assert A[[1,7]][1] eq "26b1";
    assert A[[1,9]][1] eq "54a2";
    assert A[[1,13]][1] eq "2890d1";
    assert A[[1,15]][1] eq "50a1";
    assert A[[3,3]][1] eq "19a1";
    assert A[[3,15]][1] eq "338d1";
    assert A[[4,4]][1] eq "46a1";
    assert A[[4,8]][1] eq "17a3";
    assert A[[4,12]][1] eq "14a3";
    assert A[[4,16]][1] eq "15a5";
    assert A[[4,20]][1] eq "66c1";
    assert A[[4,24]][1] eq "30a1";
    assert A[[4,32]][1] eq "210e1";
    assert A[[5,5]][1] eq "11a1";
    assert A[[8,8]][1] eq "17a1";
    assert A[[8,16]][1] eq "15a2";
    assert A[[8,24]][1] eq "30a2";
    assert A[[8,32]][1] eq "210e2";
    assert A[[12,12]][1] eq "14a1";
    assert A[[12,24]][1] eq "256a1";
    assert A[[16,16]][1] eq "15a1";
end procedure;


procedure VerifyTable2()
//compute GenD4Torsion(H) for every group in the RZB database at the appropiate level.
//We start at level 16 and then left the necessary groups up to level 32 and 64 to save time
	n := 16;
	print "Searching the level 16 torsion structures of the groups in the RZB database.";
	Tors:=AssociativeArray();
	for i in [1..#newsublist] do
    //Since the groups in RZB act on the left, we need to take their transpose.
	G := MatrixGroupTranspose(ChangeLevel(newsublist[i][3],n));
	T := TwoPowerLevelStructure(G,n);
	if T eq [] then
	T := [1,1];
	end if;
	print i, T;
	if IsDefined(Tors,T) then
    	Tors[T]:=Append(Tors[T],i);
    	else
	    Tors[T]:=[i];
	end if;
	end for;
	Lifts := Tors[[4,16]];
	for i in Tors[[8,16]] do
	Append(~Lifts,i);
	end for;
	for i in Tors[[16,16]] do
	Append(~Lifts,i);
	end for;
	Lifts := Sort(Lifts);
	Remove(~Tors,[4,16]);
	Remove(~Tors,[8,16]);
	Remove(~Tors,[16,16]);
	n := 32;
	print "Searching the level 32 torsion structures of the groups in the RZB database.";
	for i in Lifts do
	G := MatrixGroupTranspose(ChangeLevel(newsublist[i][3],n));
	T := TwoPowerLevelStructure(G,n);
	if T eq [] then
	T := [1,1];
	end if;
	i, T;
	if IsDefined(Tors,T) then
    	Tors[T]:=Append(Tors[T],i);
	    else
    	Tors[T]:=[i];
	end if;
	end for;
//Here we check which groups give curves with a point of order order 32 and check whose corresponding curves actually have points. This way we only jhave to life the curves that acttually have corresponding rational elliptic curves.
	for k in Keys(Tors) do
	if k[2] eq 32 then
	assert k in [ [ 4, 32 ],[8,32],[ 16, 32 ] ];
	end if;
	end for;
	assert Tors[[4,32]] eq [ 235, 471, 481, 496, 499, 515, 516, 517, 518, 523, 524, 534, 535, 722, 723, 724, 725 ]; //Only 235 has nontrivial rational points
	assert Tors[[ 8, 32 ]] eq [ 193, 456, 470, 475, 479, 507, 508, 509, 510 ];//Only 193 has nontrivial rational points
	assert Tors[[ 16, 32 ]] eq [ 340, 465, 467, 474, 552, 553, 554, 559, 639, 642, 644, 648 ];//None of these curves have nontrivial rational points
	print "Searching the level 64 torsion structures of the groups in the RZB database to show there are no curves with points that have a point of order 64.";
    
    n := 64;
	for i in [193,235] do
	G := MatrixGroupTranspose(ChangeLevel(newsublist[i][3],n));
	T := TwoPowerLevelStructure(G,n);
    assert T[2] eq 32;	
    end for;
    
    //Now we check that the torsion strucrues listed in Table 2 are a completel list.
    //This completes the proof of Prop 5.14.
    for k in Keys(Tors) do
             assert k in [[ 1, 1 ], [ 4, 4 ], [ 4, 8 ], [ 4, 16 ], [ 4, 32 ], [ 8, 8 ], [ 8, 16 ], [ 8, 32 ], [ 16, 16 ], [ 16, 32 ]];
	end for;
    //Each curve is now groups by its tosion structure and we check by hand using the RZB 
    //database which ones are minimal and which ones don't have rational points.
	assert Tors[[ 1, 1 ]] eq [ 1, 2, 3, 4, 5, 7, 20, 21, 22, 55, 56, 57, 177, 178, 179, 180, 439, 440, 441, 718, 719 ]; //1 is maximal
	assert Tors[[ 4, 4 ]] eq [ 6, 9, 10, 14, 15, 16, 17, 18, 19, 37, 42, 108, 123, 238, 239];//6 is maximal
	assert Tors[[ 4, 8 ]] eq [ 12, 13, 30, 31, 32, 33, 34, 40, 44, 48, 51, 52, 72, 79, 88, 89, 93, 103, 104, 115, 116, 126, 145, 157, 158, 162, 163, 164, 167, 169, 172, 173, 174, 175, 176, 226, 241, 242, 265, 266, 292, 293, 294, 295, 296, 297, 298, 299, 333, 334, 343, 368, 369, 407, 408, 419, 420, 421, 422, 434, 435, 436, 437, 548, 549, 550, 551, 555, 556, 557, 558, 563, 564, 565, 566, 571, 572, 573, 574, 616, 623, 631, 643, 645, 646, 647, 652, 653, 655, 684, 685, 704, 707, 708, 709, 710, 711, 715, 717 ]; //12, 13 are maximal
	assert Tors[[ 4, 16 ]] eq [ 36, 75, 78, 84, 85, 86, 102, 117, 118, 119, 120, 121, 122, 159, 168, 197, 199, 201, 202, 211, 212, 213, 217, 219, 223, 225, 227, 228, 229, 230, 233, 234, 236, 240, 243, 329, 330, 331, 332, 335, 336, 337, 338, 339, 341, 342, 344, 345, 346, 347, 353, 354, 355, 356, 405, 406, 427, 428, 429, 430, 431, 432, 433, 438, 466, 472, 473, 477, 488, 489, 490, 491, 492, 493, 494, 495, 497, 498, 511, 512, 513, 514, 522, 528, 530, 532, 560, 561, 562, 567, 617, 620, 621, 622, 656, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 696, 697, 698, 699, 700, 701, 702, 703, 720, 721, 726, 727 ]; //36 is maximal
	assert Tors[[ 4, 32 ]] eq [ 235, 471, 481, 496, 499, 515, 516, 517, 518, 523, 524, 534,535, 722, 723, 724, 725 ]; //235 is maximal
	assert Tors[[ 8, 8 ]] eq [ 8, 11, 23, 24, 26, 27, 28, 29, 35, 38, 39, 41, 43, 45, 46, 47, 49, 50, 53, 54, 59, 60, 61, 63, 64, 65, 66, 67, 68, 69, 70, 71, 73, 74, 76, 77, 80, 81, 82, 83, 90, 91, 94, 95, 97, 105, 106, 107, 109, 110, 111, 112, 113, 114, 124, 125, 127, 128, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 143, 144, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 160, 161, 165, 166, 170, 171, 191, 196, 209, 210, 214, 216, 218, 220, 221, 222, 224, 231, 232, 237, 250, 252, 253, 254, 256, 257, 258, 260, 261, 262, 263, 264, 267, 275, 276, 277, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 300, 301, 302, 303, 304, 307, 308, 309, 310, 311, 312, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 348, 349, 350, 351, 352, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 370, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 411, 412, 413, 414, 415, 416, 417, 418, 425, 426, 443, 461, 480, 484, 520, 525, 526, 533, 537, 538, 539, 542, 543, 545, 546, 547, 568, 569, 570, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 587, 588, 589, 590, 591, 594, 595, 596, 597, 598, 599, 600, 601, 615, 618, 619, 626, 628, 630, 632, 633, 634, 635, 637, 640, 641, 649, 650, 651, 654, 657, 658, 672, 673, 674, 675, 676, 677, 678, 679, 680, 683, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 705, 706, 712, 713, 714, 716 ]; //8, 11 are maximal
	assert Tors[[ 8, 16 ]] eq [ 25, 62, 87, 92, 96, 98, 99, 100, 101, 129, 142, 181, 182, 184, 185, 188, 190, 192, 194, 195, 198, 200, 203, 204, 205, 206, 207, 208, 215, 244, 245, 246, 247, 259, 268, 269, 270, 271, 272, 273, 274, 278, 279, 305, 306, 313, 314, 315, 316, 371, 386, 404, 409, 410, 423, 424, 442, 444, 445, 448, 451, 453, 455, 462, 463, 464, 468, 469, 476, 478, 482, 483, 485, 486, 487, 504, 505, 506, 519, 521, 527, 529, 531, 544, 592, 593, 603, 604, 605, 606, 607, 608, 609,  610, 611, 612, 613, 624, 625, 627, 629, 636, 638, 671, 681, 682 ]; //25 and 92 are maximal
	assert Tors[[ 8, 32 ]] eq [ 193, 456, 470, 475, 479, 507, 508, 509, 510 ]; //193 is maximal
	assert Tors[[ 16, 16 ]] eq [ 58, 183, 186, 187, 189, 248, 249, 251, 255, 357, 446, 447, 449, 450, 452, 454, 457, 458, 459, 460, 500, 501, 502, 503, 536, 540, 541, 586,  602, 614, 670 ]; //58 is maximal
	assert Tors[[ 16, 32 ]] eq [ 340, 465, 467, 474, 552, 553, 554, 559, 639, 642, 644, 648 ]; //None of these curves have points nontrivial rational points
end procedure;


procedure VerifyTable3()
	for k in Keys(CMTorsionArray) do
    T:= D4inftyTorsion(EllipticCurve(k):UseTmax:=false,UseRZBData:=false,UsepPrimMaps:=false, UseCMData:=false);
    assert CMTorsionArray[k] eq T;
    k,T;
    end for;
end procedure;

procedure VerifyEx31()
	D42 := DirectProduct(DihedralGroup(4),DihedralGroup(4));
    G := sub<D42| (2,4)(5,6,7,8),(1,2,3,4),(1,3)(2,4), (5,7)(6,8)>;
    H :=sub<D42| (1,3)(2,4)(5,7)(6,8) >;
    assert IsIsomorphic(SmallGroup(8,4),G/H);//Q8 is SmallGroup(8,4)
end procedure;

procedure VerifyLem32()
    G2 := SmallGroup(32,2);
    G := DirectProduct([DihedralGroup(4),DihedralGroup(4),DihedralGroup(4)]);
    t1 := G.1; t2 := G.3; t3 := G.5;
    s1 := G.2; s2 := G.4; s3 := G.6;
    H := sub<G|t1*t2,s1*t3>;
    assert IsIsomorphic(G2,H);
end procedure;

procedure VerifyLem55()
	F<s> := FunctionField(Rationals());
    Es := EllipticCurve([0,s]);
    E2s3 := EllipticCurve([0,2*s^3]);
    Es3 := EllipticCurve([0,s^3]);
    assert #TwoTorsionSubgroup(Es) eq 1 and #TwoTorsionSubgroup(E2s3) eq 1 and #TwoTorsionSubgroup(Es3) eq 2;
    assert #Factorization(DivisionPolynomial(Es,3)) eq 2 and  #Factorization(DivisionPolynomial(Es3,3)) eq 2 and #Factorization(DivisionPolynomial(E2s3,3)) eq 3;
    E1 := EllipticCurve("27a1");
    E2 := EllipticCurve("36a1");
    E3 := EllipticCurve("108a1");
    assert D4inftyTorsion(E1) eq [ 3, 3] and  D4inftyTorsion(E2) eq [ 8, 24] and  D4inftyTorsion(E3) eq [ 1, 3];
end procedure;

procedure VerifyProp58()
    G25 := GL(2,Integers(25));
    G5, pi5 := ChangeRing(G25,Integers(5));
    S125 := MaximalImagesWithGenD4Torsion([1,25]);
    S525 := MaximalImagesWithGenD4Torsion([5,25]);
    H := sub<G25| [7,0,0,1],[1,0,0,2],[1,1,0,1]>;
    for G in S125 do
    assert IsConjugateToSubgroup(G25,G,H);
    end for;
    for G in S525 do
    assert IsConjugateToSubgroup(G25,G,H);
    end for;
end procedure;

procedure VerifyProp510()
	G9 := GL(2,Integers(9));
    G3, pi3 := ChangeRing(G9,Integers(3));
    Ns3 := Normalizer(G3,GL2SplitCartan(3));
    H := sub<G9 | [1,3,0,1] , [1,0,0,2], [8,0,0,8]>;
    S33 := MaximalImagesWithGenD4Torsion([3,3]);
    S39 := MaximalImagesWithGenD4Torsion([3,9]);
    assert #S33 eq 1;
    assert IsConjugate(G3,S33[1], Ns3);
    assert IsConjugate(G9,S39[1], H);
end procedure;
             
procedure VerifyEx512()
    E := EllipticCurve("338e1");
    assert D4inftyTorsion(E) eq [3,15];
    //Therefore E has a 5-isogeny and we check that E doesn't have a 15-isogeny
    assert jInvariant(E) notin [-25/2, -349938025/8, -121945/32 , 46969655/32768 ];
end procedure;


procedure VerifyProp513()
	G4 := GL(2,Integers(4));
    G2,pi2 := ChangeRing(G4,Integers(2));
    S44 := MaximalImagesWithGenD4Torsion([4,4]);
    assert #S44 eq 1;
    assert S44[1] eq sub<G2| [1,1,0,1]>;
end procedure;

//Proposition 5.14 is verifyed in the process of verifying Table 2 so we do not create a separate procedure.


procedure VerifyProp62()
	j33 := pPrimMaps[[3,3]][1];
	j13 := pPrimMaps[[1,13]][1];
    //Construct X as the fiber product of j33 and j13
    //In the code the curve X will be denoetd by C
    R<x,y>:=PolynomialRing(Rationals(),2);
	g:=Numerator(Evaluate(j33,x)-Evaluate(j13,y));
    assert IsIrreducible(g);
    C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
    phi := iso<C -> C | [-3/X*Z^2,Y,Z], [-3/X*Z^2,Y,Z]>;
    psi := iso<C -> C | [X,-Z^2/(Y-Z),Z],[-X*Y^2*Z + X*Y*Z^2, -Y^2*Z^2 + 2*Y*Z^3 - Z^4, -Y^2*Z^2 + Y*Z^3]>;
    G := AutomorphismGroup(C,[phi,psi]);
    QC, pi := CurveQuotient(G);
	assert Genus(QC) eq 2;
    H , mp := SimplifiedModel(QC);
    J := Jacobian(H);
    assert RankBounds(J) eq 0;
    //Use the method of Chabauty and Coleman to find all the point on H.
    Pts := Chabauty0(J);
    PHI := pi*mp;
    SingPts := SingularPoints(C);
    //Pull all the points on H back to X and see that their are only singular points on X
    for P in Pts do
	assert RationalPoints(P@@PHI) subset SingPts;
	end for;
end procedure;


procedure VerifyProp63()
	j33 := pPrimMaps[[3,3]][1];
	j7 := pPrimMaps[[1,7]][1];
    //Construct X as the fiber product of j33 and j13
    //In the code the curve X will be denoetd by C and reserve X for the coodrainates
    R<x,y>:=PolynomialRing(Rationals(),2);
	g:=Numerator(Evaluate(j33,x)-Evaluate(j7,y));
    assert IsIrreducible(g);
    C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
    phi := iso<C -> C | [-3/X*Z^2,Y,Z], [-3/X*Z^2,Y,Z]>;
    psi := iso<C -> C | [X,-Z^2/(Y-Z),Z],[-X*Y^2*Z + X*Y*Z^2, -Y^2*Z^2 + 2*Y*Z^3 - Z^4, -Y^2*Z^2 + Y*Z^3]>;
    G := AutomorphismGroup(C,[phi,psi]);
    QC, pi := CurveQuotient(G);
	assert Genus(QC) eq 2;
    H , mp := SimplifiedModel(QC);
    J := Jacobian(H);
    assert RankBounds(J) eq 0;
    //Use the method of Chabauty and Coleman to find all the point on H.
    Pts := Chabauty0(J);
    PHI := pi*mp;
    SingPts := SingularPoints(C);
    //Pull all the points on H back to X and see that their are only singular points on X
    for P in Pts do
	assert RationalPoints(P@@PHI) subset SingPts;
	end for;
    //Check that elliptic curves with a 14- or 21- isogeny don't have a point of order 7 defined over Q(D4^infty)
    //These j-invariants are taken from [Table 4, 19]
    jInv := [-3^3*5^3, 3^3*5^3*17^3, -3^2*5^6/2^3, 3^3*5^3/2, -3^2*5^3*101^3/2^21, -3^3*5^3*383^3/2^7 ];
    for j in jInv do
    	E := MinimalTwist(EllipticCurveFromjInvariant(j));
        assert not HasD4inftyTorsionPoint(E,7,1: UseTmax:=false,UsepPrimMaps:=false);
     end for;
end procedure;



procedure VerifyProp64()
	j33 := pPrimMaps[[3,3]][1];
	j55 := pPrimMaps[[5,5]][1];
	//Construct X as the fiber product of j33 and j13
    //In the code the curve X will be denoetd by C and reserve X for the coodrainates
    R<x,y>:=PolynomialRing(Rationals(),2);
	g:=Numerator(Evaluate(j33,x)-Evaluate(j55,y));
	assert IsIrreducible(g);
	C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
	psi :=  iso<C -> C | [X,5*Z^2/Y,Z],[X,5*Z^2/Y,Z]>;
	G := AutomorphismGroup(C,[psi]);
	//Quotient C by the automorphism psi and construct QC
    QC,pi := CurveQuotient(G);
    //Compute the automorphism of the new curve QC call it G1
	G1 := AutomorphismGroup(QC);
    //Quotient out by the second generator of G1 to construct QC2
	QC2, pi2 := CurveQuotient(AutomorphismGroup(QC,[G1.2]));
	assert Genus(QC2) eq 1;
    //Take a singular point on C and map to to be a point on QC2
    P2 := C ! [0,1,0];
    R2 := pi(P2);
	Q22 := pi2(R2);
    //Use the point Q22 to get a Weierstrass model for QC2 Call it E1
	E, pi3 := EllipticCurve(QC2,Q22);
	//Show E is a quadratic twist of 15a8 by a rational square so they are isomorphic over Q
    E1 := EllipticCurve("15a8");
	TF, d := IsQuadraticTwist(E1,E);
	assert IsSquare(d);
    //Thus curve E has M--W Group isomorphic to Z/4Z we construct a point of order 4. Call it P.
	a:= Roots(DivisionPolynomial(E,4) div DivisionPolynomial(E,2))[2][1];
	A := aInvariants(E);
	_<y> := PolynomialRing(Rationals());
	g := y^2 + A[1]*a*y + A[3]*y - (a^3 + A[2]*a^2 + A[4]*a + A[5]);
	b := Roots(g)[1][1];
	P := E![a,b];
	assert P ne Identity(E) and 2*P ne Identity(E) and 4*P eq Identity(E);
    //Create a list of the points on E and pull them back to QC2
	PtsE := [n*P : n in [0..3]];
	InvQC2 := [RationalPoints((Q)@@pi3): Q in PtsE];
	PtsQC2 := [];
	for list in InvQC2 do
	for Q in list do
	if Q notin PtsQC2 then
	Append(~PtsQC2,Q);
	end if;
	end for;
	end for;
    //Pull the points on QC2 back to QC
    InvQC := [RationalPoints((Q)@@pi2): Q in PtsQC2];
    PtsQC := [];
	for list in InvQC do
	for Q in list do
	if Q notin PtsQC then
	Append(~PtsQC,Q);
	end if;
	end for;
	end for;
    //Pull the points on QC back to C
	InvC := [RationalPoints((Q)@@pi): Q in PtsQC];
	PtsC := [];
	for list in InvC do
	for Q in list do
	if Q notin PtsC then
	Append(~PtsC,Q);
	end if;
	end for;
	end for;
    assert PtsC subset SingularPoints(C);
end procedure;


procedure VerifyProp65()
	j33 := pPrimMaps[[3,3]][1];
	j5 := pPrimMaps[[1,5]][1];
    //Construct the curve X as the fiber product of j33 and j5
    //Use C to denote the curve instead of X
    R<x,y>:=PolynomialRing(Rationals(),2);
	g:=Numerator(Evaluate(j33,x)-Evaluate(j5,y));
	assert IsIrreducible(g);
	C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
    assert Genus(C) eq 1;
    SingPtsC := SingularPoints(C);
    P := C![0,1,0];
    //Use a singular point to get a Weierstrass equation for C and call the new curve E
    E, pi := EllipticCurve(C,P);
	//Show the rank of E(Q) is 0
    assert Rank(E) eq 0;
    //Take every points on E and pull it back to C and then plug the nonsingular points in to either
    //j-map to see that there are only 2 Qbar-isomorphism classes of curves with torsion structure Z/3Z x Z/15Z
    T , map := TorsionSubgroup(E);
    for P in T do
	for Q in [R : R in RationalPoints(map(P) @@ pi)] do
	if not Q in SingPtsC then
	assert Evaluate(j33,Q[1]/Q[3]) eq Evaluate(j5,Q[2]/Q[3]);
    assert Evaluate(j33,Q[1]/Q[3]) in [1331/8,-1680914269/32768];
	end if;
	end for;
	end for;
end procedure;

procedure VerifyProp66()
	//We first rule out possibility (1) for Z/4Z x Z/8Z
    F<t> := FunctionField(Rationals());
	//Construct the generic elliptic curve with a raitonal 10-isogeny. The j-invariant is taken from Table 3 of [18]
    E := EllipticCurveFromjInvariant((t^6-4*t^5+16*t+16)^3*(t+1)^-2*(t-4)^-1*t^-5);
    G, map := TwoTorsionSubgroup(E);
	P := map(G.1);
	//Construct E2, the curve that is 2-isogenous to E.
    E2 := Image(TwoIsogeny(P));
	//Computing the ratio of their disriminants R, we see that (1) is satifies if and only R is rational a square.
    D := Discriminant(E);
	D2 := Discriminant(E2);
	R := D/D2;
    assert R eq (t+1)/(t^6-4*t^5);
    //Define the curve S^2 = R
	A<T,S> := AffineSpace(Rationals(),2);
	CS := Curve(A, (T+1) - (T^6-4*T^5)*S^2);
    assert Genus(CS) eq 1;
	Pt := CS![-1,0];
	Er, mapEr := EllipticCurve(ProjectiveClosure(CS),Pt);
	assert Rank(Er) eq 0;
	T2,mapT2 := TwoTorsionSubgroup(Er);
	Pts := [mapT2(Q) : Q in T2];
	for P in Pts do
	assert RationalPoints(P@@mapEr) subset {@ ProjectiveClosure(CS)! [-1 , 0 , 1],  ProjectiveClosure(CS)![0 , 1 , 0],  ProjectiveClosure(CS)![1 , 0 , 0] @};
	end for;
    //Next we rule out possibility (1) for Z/8Z x Z/8Z
    A<T,S> := AffineSpace(Rationals(),2);
	//The curve is the same elliptic curve but it is twited by -1.
    CS1 := Curve(A, (T+1) + (T^6-4*T^5)*S^2);
    assert Genus(CS) eq 1;
	Pt := CS![-1,0];
	Er, mapEr := EllipticCurve(ProjectiveClosure(CS),Pt);
	assert Rank(Er) eq 0;
	T2,mapT2 := TwoTorsionSubgroup(Er);
	Pts := [mapT2(Q) : Q in T2];
	for P in Pts do
	assert RationalPoints(P@@mapEr) subset {@ ProjectiveClosure(CS)! [-1 , 0 , 1],  ProjectiveClosure(CS)![0 , 1 , 0],  ProjectiveClosure(CS)![1 , 0 , 0] @};
	end for;
end procedure;

procedure VerifyProp67()
	J := [1331/8,-1680914269/32768];
    for j in J do
    E := MinimalTwist(EllipticCurveFromjInvariant(j));
    assert D4inftyTorsion(E:UseTmax:=false,UseRZBData:=false,UsepPrimMaps:=false, UseCMData:=false) eq [3,15];
    end for;
end procedure;



procedure VerifyProp610()
	//First we show that there are no elliptic curves with full 3- and full 8- torsion defined over Q(D4^ifnty)
    j881 := pPrimMaps[[8,8]][1];
	j882 := pPrimMaps[[8,8]][2];
	j33 := pPrimMaps[[3,3]][1];
    //We construct the two corresponding modular curves as the fiber product of j33 with the two relevant j-maps
    R<x,y>:=PolynomialRing(Rationals(),2);
	g1:=Numerator(Evaluate(j33,x)-Evaluate(j881,y));
	g2:=Numerator(Evaluate(j33,x)-Evaluate(j882,y));
	assert IsIrreducible(g1);
    assert IsIrreducible(g2);
    C1<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g1));
	C2<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g2));
    assert Genus(C1) eq 1;
    assert Genus(C2) eq 1;
    //Use the singule point [0,1,0] on both curves to get a Weierstrass equesiton for
    //C1 and C2. Call them E1 and E2 respectively
    P1 := C1 ![0,1,0];
    P2 := C2 ![0,1,0];
	E1, pi1 := EllipticCurve(C1,P1); //Cremona:36a2
	E2, pi2 := EllipticCurve(C2,P2); //Cremona:48a1
    assert Rank(E1) eq 0;
    assert Rank(E2) eq 0;
    //Since both curves have rank 0, all their points are torsion points
    //Construct to torsion subgroups and pull them back to C1 and C2 respectively
	T1, map1 := TorsionSubgroup(E1);
	T2, map2 := TorsionSubgroup(E2);
    //Show all the rational points on C1 are singular
    for Q in T1 do
	assert RationalPoints(map1(Q) @@ pi1) subset SingularPoints(C1);
	end for;
    //Show all the rational points on C2 are singular
    for Q in T2 do
	assert RationalPoints(map2(Q) @@ pi2) subset SingularPoints(C2);
	end for;
    //Next we show that there are no elliptic curves with full 3- and a point of order 16 defin over Q(D4^infty)
    j8 := pPrimMaps[[4,16]][1]; //the modular curve X36 in RZB is the classical modular curve X0(8)
    g:=Numerator(Evaluate(j33,x)-Evaluate(j8,y));
    assert IsIrreducible(g);
    C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
    assert Genus(C) eq 3;
    assert IsHyperelliptic(C);
    phi := iso<C -> C | [-3/X*Z^2,Y,Z], [-3/X*Z^2,Y,Z]>;
    G := AutomorphismGroup(C,[phi]);
    H, pi := CurveQuotient(G);
    J := Jacobian(H);
    assert RankBounds(J) eq 0;
    PtsH := Chabauty0(J);
    for Q in PtsH do
    assert RationalPoints(Q@@pi) subset {@ C![0 , -4 , 1], C![0 , 0 , 1], C![0 , 4 , 1], C![0 , 1 , 0], C![1 , 0 , 0] @};
    end for;
end procedure;


// List of the rational functions that appear in Table 4 (some of which are constant functions)
F<t>:=FunctionField(Rationals());
Tparam := AssociativeArray();
Tparam[[ 1, 1 ]] := <t>;
Tparam[[ 1, 3 ]] := <(27*t^4 + 756*t^3 + 7290*t^2 + 26244*t + 19683)/t^3>;
Tparam[[ 1, 5 ]] := <25 * (t^2+10*t+5)^3/t^5>;
Tparam[[ 1, 7 ]] := <(t^2-t+1)^3*(t^6-11*t^5+30*t^4-15*t^3-10*t^2+5*t+1)^3*(t-1)^-7*t^-7*(t^3-8*t^2+5*t+1)^-1 >;
Tparam[[ 1, 9 ]] := <(t^36 - 36*t^35 + 594*t^34 - 6000*t^33 + 41859*t^32 - 215892*t^31 + 860466*t^30 - 2733948*t^29 + 7083369*t^28 - 15215576*t^27 + 27417258*t^26 - 41737032*t^25 + 53766648*t^24 - 58258332*t^23 + 52042986*t^22 - 36346716*t^21 + 16629516*t^20 + 241380*t^19 - 9621660*t^18 + 11057688*t^17 - 7521093*t^16 + 2914848*t^15 + 178308*t^14 - 1208520*t^13 + 948843*t^12 - 375948*t^11 + 15714*t^10 + 78304*t^9 - 48717*t^8 + 11124*t^7 + 2262*t^6 - 2124*t^5 + 378*t^4 + 84*t^3 - 36*t^2 + 1)/(t^27 - 18*t^26 + 144*t^25 - 702*t^24 + 2376*t^23 - 5985*t^22 + 11688*t^21 - 18135*t^20 + 22662*t^19 - 22923*t^18 + 18711*t^17 - 12168*t^16 + 6132*t^15 - 2259*t^14 + 522*t^13 - 27*t^12 - 27*t^11 + 9*t^10 - t^9)>;
Tparam[[ 1, 13 ]] := <(t^2 - t + 1)^3*( t^12 - 9*t^11 + 29*t^10 - 40*t^9 + 22*t^8 - 16*t^7 + 40*t^6 - 22*t^5 - 23*t^4 + 25*t^3 - 4*t^2 - 3*t + 1)^3 * (t - 1)^(-13)*t^(-13)*(t^3 - 4*t^2 + t +1)^-1>;
Tparam[[ 1, 15 ]] := <F!-25/2, F!-349938025/8, F!-121945/32 , F!46969655/32768 >;
Tparam[[ 3, 3 ]] := <(27*(t+1)^3*(t-3)^3)/t^3>;
Tparam[[ 3, 15 ]] := <F!1331/8,F!-1680914269/32768>;
Tparam[[ 4, 4 ]] := <t^3/(t + 16)>;
Tparam[[ 4, 8 ]] := <(t^6 + 48*t^4 + 768*t^2 + 4096)/t^2, (t^6 - 144*t^4 + 6912*t^2 - 110592)/(t^2 - 64)>;
Tparam[[ 4, 12 ]] := <(t+6)^3*(t^3+18*t^2+84*t+24)^3*t^-1*(t+8)^-3*(t+9)^-2>;
Tparam[[ 4, 16 ]] := <(t^12 - 48*t^10 + 816*t^8 - 5632*t^6 + 13056*t^4 - 12288*t^2 + 4096)/(t^4 - 16*t^2)>;
Tparam[[ 4, 20 ]] := <(t^6-4*t^5+16*t+16)^3/((t+1)^2*(t-4)*t^5)>;
Tparam[[ 4, 24 ]] := <(t^2-3)^3*(t^6-9*t^4+3*t^2-3)^3/(t^4*(t^2-9)*(t^2-1)^3)>;
Tparam[[ 4, 32 ]] := <(t^48 - 24*t^46 + 228*t^44 - 1064*t^42 + 2322*t^40 - 840*t^38 - 5676*t^36 + 8712*t^34 + 1455*t^32 - 16240*t^30 + 14664*t^28 + 9456*t^26 - 21892*t^24 + 9456*t^22 + 14664*t^20 - 16240*t^18 + 1455*t^16 + 8712*t^14 - 5676*t^12 - 840*t^10 + 2322*t^8 - 1064*t^6 + 228*t^4 - 24*t^2 + 1)/(t^32 - 8*t^30 + 12*t^28 + 8*t^26 - 26*t^24 + 8*t^22 + 12*t^20 - 8*t^18 + t^16)>;
Tparam[[ 5, 5 ]] := <(t^30 + 30*t^29 + 465*t^28 + 4900*t^27 + 39150*t^26 + 252000*t^25 + 1360000*t^24 + 6337500*t^23 + 26090625*t^22 + 96593750*t^21 + 325846875*t^20 + 1010625000*t^19 + 2898984375*t^18 + 7721718750*t^17 + 19152890625*t^16 + 44323437500*t^15 + 95764453125*t^14 + 193042968750*t^13 + 362373046875*t^12 + 631640625000*t^11 + 1018271484375*t^10 + 1509277343750*t^9 + 2038330078125*t^8 + 2475585937500*t^7 + 2656250000000*t^6 + 2460937500000*t^5 + 1911621093750*t^4 + 1196289062500*t^3 + 567626953125*t^2 + 183105468750*t + 30517578125)/(t^25 + 25*t^24 + 325*t^23 + 2875*t^22 + 19250*t^21 + 103125*t^20 + 456875*t^19 + 1709375*t^18 + 5475000*t^17 + 15140625*t^16 + 36321875*t^15 + 75703125*t^14 + 136875000*t^13 + 213671875*t^12 + 285546875*t^11 + 322265625*t^10 + 300781250*t^9 + 224609375*t^8 + 126953125*t^7 + 48828125*t^6 + 9765625*t^5)>;
Tparam[[ 8, 8 ]] := <(64*t^6 + 576*t^4 + 1728*t^2 + 1728)/(t^4 - 2*t^2 + 1), (t^6 - 48*t^4 + 768*t^2 - 4096)/t^2>;
Tparam[[ 8, 16 ]] := <(256*t^12 - 768*t^10 + 1536*t^8 - 1792*t^6 + 1536*t^4 - 768*t^2 + 256)/(t^8 - 2*t^6 + t^4), (-4*t^24 + 720*t^22 - 44808*t^20 + 1057680*t^18 - 6090684*t^16 + 15714720*t^14 - 21291632*t^12 + 15714720*t^10 - 6090684*t^8 + 1057680*t^6 - 44808*t^4 + 720*t^2 - 4)/(t^22 + 6*t^20 + 13*t^18 + 8*t^16 - 14*t^14 - 28*t^12 - 14*t^10 + 8*t^8 + 13*t^6 + 6*t^4 + t^2)>;
Tparam[[ 8, 24 ]] := <(1/64*t^24 - 9/8*t^23 + 585/16*t^22 - 5643/8*t^21 + 284841/32*t^20 - 605259/8*t^19 + 6718797/16*t^18 - 10181457/8*t^17 - 11997153/64*t^16 + 64644075/4*t^15 - 153327303/8*t^14 - 1444643343/4*t^13 + 31741088463/16*t^12 - 9567551763/4*t^11 - 108426844863/8*t^10 + 175634578935/4*t^9 + 5659193923407/64*t^8 - 5680109894877/8*t^7 + 15711151677957/16*t^6 + 18233184917241/8*t^5 - 225428092751799/32*t^4 - 55488939877503/8*t^3 + 856650746267505/16*t^2 - 676525954008429/8*t + 2991025682883441/64)/(t^18 - 54*t^17 + 1281*t^16 - 17424*t^15 + 147444*t^14 - 767592*t^13 + 2004388*t^12 + 2006928*t^11 - 37780626*t^10 + 139766076*t^9 - 205242594*t^8 - 233825904*t^7 + 1766689156*t^6 - 4060927080*t^5 + 5446948500*t^4 - 4640490000*t^3 + 2479865625*t^2 - 760893750*t + 102515625)>;
Tparam[[ 8, 32 ]] := <(256*t^48 - 1536*t^46 + 3648*t^44 - 4256*t^42 + 3042*t^40 - 3090*t^38 + 14421/4*t^36 - 14751/8*t^34 + 171375/256*t^32 - 70855/64*t^30 + 235833/512*t^28 + 79071/1024*t^26 + 3134807/16384*t^24 + 79071/16384*t^22 + 235833/131072*t^20 - 70855/262144*t^18 + 171375/16777216*t^16 - 14751/8388608*t^14 + 14421/67108864*t^12 - 1545/134217728*t^10 + 1521/2147483648*t^8 - 133/2147483648*t^6 + 57/17179869184*t^4 - 3/34359738368*t^2 + 1/1099511627776)/(t^40 - 4*t^38 + 11/2*t^36 - 11/4*t^34 - 9/64*t^32 + 39/64*t^30 - 83/512*t^28 - 27/1024*t^26 + 611/32768*t^24 - 27/16384*t^22 - 83/131072*t^20 + 39/262144*t^18 - 9/4194304*t^16 - 11/4194304*t^14 + 11/33554432*t^12 - 1/67108864*t^10 + 1/4294967296*t^8)>;
Tparam[[ 12, 12 ]] := <(1/729) * (t - 11)^-3 * (t - 8)^-6 * (t + 1)^-3 * (t^3 - 15*t^2 - 33*t + 955)^3 * (t^3 - 15*t^2 + 75*t - 233)^3>;
Tparam[[ 12, 24 ]] := <F!8000>;
Tparam[[ 16, 16 ]] := <(16*t^24 + 672*t^20 + 9456*t^16 + 45248*t^12 + 9456*t^8 + 672*t^4 + 16)/(t^20 - 4*t^16 + 6*t^12 - 4*t^8 + t^4)>;


procedure VerifyTable4()
	//There are only 5 cases that are left to be verifyied they are (Z/3Z x Z/15Z), (Z/4Z x Z/24Z), (Z/8Z x Z/24Z), (Z/12Z x Z/1Z), and (Z/12Z x Z/24Z)
	
    //The case of (Z/3Z x Z/15Z) was already verified in Proposition 6.5 and 6.7.
    
    //The first part of case of (Z/4Z x Z/24Z) (the second part is cited in the literature).
    F<t> := FunctionField(Rationals());
    //Let E be the generic elliptic curve with a 6-isogney. The j-map is taken from Table 3 [18]
    E := EllipticCurveFromjInvariant((t+6)^3*(t^3+18*t^2+84*t+24)^3/(t*(t+8)^3*(t+9)^2));
    G, map := TwoTorsionSubgroup(E);
	P := map(G.1);
	//Construct E2, the curve that is 2-isogenous to E.
    E2 := Image(TwoIsogeny(P));
	//Computing the ratio of their disriminants R.
    D := Discriminant(E);
	D2 := Discriminant(E2);
	R := D/D2;
	assert R eq (t + 9)/(t^4 + 24*t^3 + 192*t^2 + 512*t);
    //Define the curve S^2 = R to show that there are no nonsingular rational values of t that make R a square.
    A<T,S> := AffineSpace(Rationals(),2);
	CS := Curve(A, (T+9) - (T^4 + 24*T^3 + 192*T^2 + 512*T)*S^2);
    assert Genus(CS) eq 1;
	Pt := CS![-9,0];
    //Find a Weierstrass mode for CS1 using Pt1
	Er, mapEr := EllipticCurve(ProjectiveClosure(CS),Pt);
	assert Rank(Er) eq 0;
	T2,mapT2 := TwoTorsionSubgroup(Er);
	Pts := [mapT2(Q) : Q in T2];
	for P in Pts do
	assert RationalPoints(P@@mapEr) subset {@ ProjectiveClosure(CS)! [-9, 0 , 1],  ProjectiveClosure(CS)![0 , 1 , 0],  ProjectiveClosure(CS)![1 , 0 , 0] @};
	end for;
    
    //The first part of case of (Z/8Z x Z/24Z) (the second part is cited in the literature).
    //Let E be the generic elliptic curve with a 6-isogney. The j-map is taken from Table 3 [18]
    E := EllipticCurveFromjInvariant((t+6)^3*(t^3+18*t^2+84*t+24)^3/(t*(t+8)^3*(t+9)^2));
    G, map := TwoTorsionSubgroup(E);
	P := map(G.1);
	//Construct E2, the curve that is 2-isogenous to E.
    E2 := Image(TwoIsogeny(P));
	//Computing the ratio of their disriminants R.
    D := Discriminant(E);
	D2 := Discriminant(E2);
	R := D/D2;
	assert R eq (t + 9)/(t^4 + 24*t^3 + 192*t^2 + 512*t);
    //Define the curve S^2 = R to show that there are no nonsingular rational values of t that make R a square.
    A<T,S> := AffineSpace(Rationals(),2);
	CS1 := Curve(A, (T+9) + (T^4 + 24*T^3 + 192*T^2 + 512*T)*S^2);
    assert Genus(CS1) eq 1;
	Pt1 := CS1![-9,0];
    //Find a Weierstrass mode for CS1 using Pt1
	Er1, mapEr1 := EllipticCurve(ProjectiveClosure(CS1),Pt1);
	assert Rank(Er1) eq 0;
	T21,mapT21 := TwoTorsionSubgroup(Er1);
	Pts1 := [mapT21(Q) : Q in T21];
	for P in Pts1 do
	assert RationalPoints(P@@mapEr1) subset {@ ProjectiveClosure(CS1)! [-9, 0 , 1],  ProjectiveClosure(CS1)![0 , 1 , 0],  ProjectiveClosure(CS1)![1 , 0 , 0] @};
	end for;

	//The case of (Z/12Z x Z/12Z)
    j44 := pPrimMaps[[4,4]][1];
	j33 := pPrimMaps[[3,3]][1];
    //Consturct the modular curve X as the fiber product of j33 and j44
	R<x,y>:=PolynomialRing(Rationals(),2);
	g:=Numerator(Evaluate(j33,x)-Evaluate(j44,y));
	C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
	assert Genus(C) eq 0;
	//Use the universal mapping property of fiber products and the parameterization of C to compute the j-map
    Co,phi1:=Conic(C);
	P1:=Curve(ProjectiveSpace(Rationals(),1));
	phi2:=Parametrization(Co,P1);
	phi3:=Expand(phi2*Inverse(phi1));
	phi3polys:=DefiningPolynomials(phi3);
    h:=Evaluate(j33,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 case of (Z/12Z x Z/24Z)
	j33 := pPrimMaps[[3,3]][1];
    j481 := pPrimMaps[[4,8]][1];
	j482 := pPrimMaps[[4,8]][2];
    //Construct the corresponding modular curves C1 and C2
	R<x,y> := PolynomialRing(Rationals(),2);
    g1 := Numerator(Evaluate(j33,x)-Evaluate(j481,y));
	g2 := Numerator(Evaluate(j33,x)-Evaluate(j482,y));
    C1<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g1));
    C2<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g2));
    P1 := C1 ! [0,1,0];
    P2 := C2 ! [0,1,0];
    E1,pi1 := EllipticCurve(C1,P1);
    E2,pi2 := EllipticCurve(C2,P2);
    assert Rank(E1) eq 0 and Rank(E2) eq 0;
    T1, map1 := TorsionSubgroup(E1);
    T2, map2 := TorsionSubgroup(E2);
	//Pullback the points in T1 to C1 and see that all the rational nonsingular points correspond to j = 8000
    for P in T1 do
	list := RationalPoints(map1(P) @@ pi1);
	for Q in list do
	if not Q in SingularPoints(C1) then
	assert Evaluate(j33, Q[1]/Q[3]) eq Evaluate(j481, Q[2]/Q[3]) and  Evaluate(j33, Q[1]/Q[3]) eq 8000;
	end if;
	end for;
	end for;
    //Pullback the points in T1 to C1 and see that all the rational nonsingular points are cusps
    for P in T2 do
	list := RationalPoints(map2(P) @@ pi2);
	for Q in list do
	if not Q in SingularPoints(C2) then
	assert Evaluate(Denominator(j33), Q[1]/Q[3]) eq 0 and Evaluate(Denominator(j482), Q[2]/Q[3]) eq 0;
	end if;
	end for;
	end for;
end procedure;


function D4inftyTorsionFast(E)
    if (Type(E) eq RngIntElt or Type(E) eq FldRatElt) then
        if E eq 0 then print "For j-invariant 0 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 HasCm(E) and jE ne 0 then
    for k in Keys(CMTorsionArray) do
        	if IsQuadraticTwist(E, EllipticCurve(k)) then
            return CMTorsionArray[k];
            end if;
        end for;
    end if;
    if jE eq 0 then
        return D4inftyTorsion(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;


// Verify that the parameterizations for torsion group T produce a curve E/Q for which E(Q(D4^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

//Unforunately it is computationally unfeasable to verify the cases of T = [4,8], [4,16], [4,32], [8,16], or [8,32] this way and in each case we
//are limited to checking that T contained either the full 4-torsion or full 8- torsion respectively. The good news is that these parameterization will
//be check for curves of small conductor (less than 500) in ParamCheck2 and there are many curves with small conductor and torsion subgroup isomorphic
//to [4,8], [4,16], [4,32], [8,16], or [8,32]. For example, there are 333, 130, 2, 72, and 1 curve respectively with onductor less then 100 and torsion
//subgroup isomorphic to [4,8], [4,16], [4,32], [8,16], or [8,32].

function ParamCheck1(n:UseTmax := false, UsepPrimMaps := false)
    printf "Sanity checking each parameterization at %o random j-invariants, this will take a while (possibly hours), be patient...\n", n;
    K := Sort([k : k in Keys(Tparam)]);
    for T in K do
    	if not T in [[4,8],[4,16],[4,32],[8,16],[8,32]] then
        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,50)/Random(1,50);
                if Evaluate(Denominator(f),r) ne 0 then
                    E:=EllipticCurveFromjInvariant(Evaluate(f,r));
                    Etor:=D4inftyTorsion(E: UseTmax := false, UsepPrimMaps := false);
                    if not TorsionContains(Etor,T) then
                        printf "Verification for T = %o failed using function f = %o with t = %o yielding E = %o, which has torsion group %0\n", T, f, r, E, Etor;
                        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 if;
//Here we check that curves that are supposed to have torsion subgroup [4,8],[4,16], or [4,32] at least contain full 8-torsion.
        if  T in [[4,8],[4,16],[4,32]] then
        for f in Tparam[T] do
            k:=0; m:=0;
            for i:=1 to n do
                r:=Random(1,50)/Random(1,50);
                if Evaluate(Denominator(f),r) ne 0 then
                    E:=EllipticCurveFromjInvariant(Evaluate(f,r));
                    if not HasD4inftyFullTorsion(E,2,2) then
                        printf "Verification for T = %o failed using function f = %o with t = %o yielding E = %o", T, f, r, E;
                        assert false;
                    end if;
                    k +:= 1;
                    if HasD4inftyFullTorsion(E,2,2)  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 curves containing [ 4, 4 ] \n", T, f, k, m;
        end for;
        end if;
//Here we check that curve that are supposed to have torsion subgroup [8,16] or [8,32] at least contain full 8-torsion.
        if  T in [[8,16],[8,32]] then
        for f in Tparam[T] do
            k:=0; m:=0;
            for i:=1 to n do
                r:=Random(1,50)/Random(1,50);
                if Evaluate(Denominator(f),r) ne 0 then
                    E:=EllipticCurveFromjInvariant(Evaluate(f,r));
                    if not HasD4inftyFullTorsion(E,2,3) 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 HasD4inftyFullTorsion(E,2,3) 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 curves containing [ 8, 8 ]\n", T, f, k, m;
        end for;
        end if;
        
    end for;
    printf "Successfully checked each parameterization at %o random points", n;
    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(D4^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:=D4inftyTorsion(E: UseTmax:=false,UseRZBData:=false,UsepPrimMaps:=false, UseCMData:=false);
            printf "%o %o\n", CremonaReference(E), T;
            assert D4inftyTorsionFast(E) eq T;
            if jInvariant(E) eq 0 then assert T in {[3,3], [8,24], [1,3]}; end if;
        end for;
    end for;
    printf "Successfully checked E(Q(D4^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 Table 1 ***";
    VerifyTable1();
    print "*** Table 1 Verified ***";

    print "*** Verifying Table 2 ***";
    VerifyTable2();
    print "*** Table 2 Verified ***";

    print "*** Verifying Table 3 ***";
    VerifyTable3();
    print "*** Table 3 Verified ***";

    print "*** Verifying Example 3.1 ***";
    VerifyEx31();
    print "*** Example 3.1 Verified ***";

    print "*** Verifying Lemma 3.2 ***";
    VerifyLem32();
    print "*** Lemma 3.2 Verified ***";

    print "*** Verifying Lemma 5.5 ***";
    VerifyLem55();
    print "*** Lemma 5.5 Verified ***";

    print "*** Verifying Proposition 5.8 ***";
    VerifyProp58();
    print "*** Proposition 5.8 Verified ***";

    print "*** Verifying Proposition 5.10 ***";
    VerifyProp510();
    print "*** Proposition 5.10 Verified ***";

    print "*** Verifying Example 5.12 ***";
    VerifyEx512();
    print "*** Example 5.12 Verified ***";

    print "*** Verifying Lemma 6.2 ***";
    VerifyProp62();
    print "*** Proposition 6.2 Verified ***";

    print "*** Verifying Remark 6.3 ***";
    VerifyProp63();
    print "*** Remark 6.3 Verified ***";

    print "*** Verifying Proposition 6.4 ***";
    VerifyProp64();
    print "*** Proposition 6.4 Verified ***";

    print "*** Verifying Proposition 6.5 ***";
    VerifyProp65();
    print "*** Proposition 6.5 Verified ***";

    print "*** Verifying Proposition 6.6 ***";
    VerifyProp66();
    print "*** Proposition 6.6 Verified ***";

    print "*** Verifying Proposition 6.10 ***";
    VerifyProp610();
    print "*** Proposition 6.10 Verified ***";

    print "*** Verifying Table 4 ***";
    VerifyTable4();
    print "*** Table 4 ***";

    print "*** Sanity checking parametrization in Table 4 ***";
    ParamCheck1(3);
    ParamCheck2(500);
    
    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 a while.";
print "The function D4inftyTorsion(E) computes E(Q(D4^infty))_tor for any elliptic curve over E/Q.";
print "The function D4inftyTorsionFast(E) is a faster version that uses the parametrizations of Table 4.";




