/*
    Magma code related to the paper "Torsion subgroups of rational elliptic curves over the compositum of all extensions of generalized $A_4$-type"
    by H. Daniels, M. Derickx, and J. Hatley
    Copyright (c) 2018 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)
	Much of this code has been taken from https://math.mit.edu/~drew/Q3tor/Q3tor.m and modified to fit our needs.


    Run VerifyAll() to verify all the computations cited in the paper (takes about 24 hours on a fast machine).

*/

//Load RZB Database of groups
//load "~/Dropbox (Amherst College)/Research/Gen-A4/RZB Data/gl2data";
//load "RZB Data/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 to 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;

// Subgroup of GL(2,Z/7Z) which is labeled 7NS.2.1 in Sutherland's Sigma paper
function NS1(x)
    R:=Integers(7);
    return sub<GL(2,R) | [[2,0,0,4]], [[0,1,4,0]]>;
end function;

// Subgroup of GL(2,Z/7Z) which is labeled 7NS.3.1 in Sutherland's Sigma paper
function NS2(x)
    R:=Integers(7);
    return sub<GL(2,R) | [[3,0,0,5]], [[0,1,4,0]]>;
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 A4-type, if it returns false then it may or may not be.
function NotGenA4Poly(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;

function LowerCentralpSeries(G,p)
	D := DerivedGroup(G);
	Gens := {@ g^p : g in G | g^p ne G.0 @} join {@ g : g in Generators(D) @};
	return sub<G| Gens>;
end function;

// Fast algorithm to check whether a number field or group is of generalized A4-type via Theorem 3.2
function IsGenA4(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 NotGenA4Poly(DefiningPolynomial(G)) then return false; else return $$(GaloisGroup(G)); end if;
    end if;
    if Type(G) eq RngUPolElt then
        if NotGenA4Poly(G) then return false; else return $$(GaloisGroup(G)); end if;
    end if;
    return #LowerCentralpSeries(LowerCentralpSeries(G,3),2) eq 1;
end function;


// Creates the group of units in Z/nZ
function UnitGroup(n)
    return GeneralLinearGroup(1,Integers(n));
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 A4-type
function GenA4Divisor(G)
    if IsGenA4(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 IsGenA4(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(A4^infty))
function GenA4Torsion(G)
    A := FixModule(GenA4Divisor(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(A_4^infy))[n]
// is isomorphic to T.  The returned groups will each be reduced modulo m=GL2Level(G) (unless m=1)
function MaximalImagesWithGenA4Torsion(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|GenA4Torsion(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;

F<t> := FunctionField(Rationals());
pPrimMap:=AssociativeArray();
pPrimMap[[ 1, 1 ]] := <t>;
pPrimMap[[ 1, 3 ]] := <(t+27)*(t+3)^3/t>;
pPrimMap[[ 1, 5 ]] := <(t^4-12*t^3+14*t^2+12*t+1)^3/(t^5*(t^2-11*t-1))>;
pPrimMap[[ 1, 7 ]] := <(t^2+13*t+49)*(t^2+5*t+1)^3/t>;
pPrimMap[[ 1, 9 ]] := <t^3*(t^3-24)^3/(t^3-27)>;
pPrimMap[[ 1, 13 ]] := <(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))>;
pPrimMap[[ 2, 2 ]] :=  <t^3/(t + 16)>;
pPrimMap[[ 2, 4 ]] := <(t^6 - 144*t^4 + 6912*t^2 - 110592)/(t^2 - 64)>;
pPrimMap[[ 2, 8 ]] := <(t^12 - 48*t^10 + 816*t^8 - 5632*t^6 + 13056*t^4 - 12288*t^2 + 4096)/(t^4 - 16*t^2)>;
pPrimMap[[ 2, 16 ]] := <(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)>;
pPrimMap[[ 3, 3 ]] := < 27*(t+1)^3*(t+3)^3*(t^2+3)^3/(t^3*(t^2+3*t+3)^3) >;
pPrimMap[[ 3, 9 ]] := < 0 >;
pPrimMap[[ 4, 4 ]] := <t^2 + 1728, (-64*t^12 + 672*t^10 - 2364*t^8 + 2828*t^6 - 591/4*t^4 + 21/8*t^2 - 1/64)/(t^10 + t^8 + 3/8*t^6 + 1/16*t^4 + 1/256*t^2) >;
pPrimMap[[ 4, 8 ]] := <(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)>;
pPrimMap[[ 4, 16 ]] := <(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)>;
pPrimMap[[ 8, 8 ]] := <(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)>;





//Given E/Q, determines if E(Q(A4^infy))[p^n]=E[p^n].
function HasA4inftyFullTorsion(E,p,n)
    assert IsPrime(p) and n ge 1;
    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 IsGenA4(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 IsGenA4(L) then return false; end if;
    end for;
    return true;
end function;

//Given E/Q, determines if E(Q(A4^infy))(p) contains a point of order p^n.

function HasA4inftyTorsionPoint(E,p,n)
    assert IsPrime(p) and n ge 1;
    if not p in [2,3,5,7,13] then return false; 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 IsGenA4(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 IsGenA4(L) then return true; end if;
    end for;
    return false;
end function;

//Given E/Q, determines  E(Q(A4^infy))(p)

function pPrimaryA4inftyTorsion(E,p)
    E:=WeierstrassModel(IntegralModel(E));
    m:=0;
    while true do
        if not HasA4inftyFullTorsion(E,p,m+1) then break; end if;
        m:=m+1;
    end while;
    n:=m;
    while true do
        if not HasA4inftyTorsionPoint(E,p,n+1) then break; end if;
        n:=n+1;
    end while;
    return [p^m,p^n];
end function;

//Given E/Q, determines E(Q(A4^infy))_tor

function A4inftyTorsion(E)
    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:=pPrimaryA4inftyTorsion(E,p);
        T[1]:=T[1]*Tp[1]; T[2]:=T[2]*Tp[2];
    end for;
    return T;
end function;

//Computes E(Q(A4))_tor for every elliptic curve with conductor less than maxN and returns an array
//where the keys are torsion strucutures and vectors contain the cremona labels of curves with that
//torsion strucutre.
function A4inftyTorsionExamples(maxN)
    examples:=AssociativeArray();
    D:=CremonaDatabase();
    for N in [11..maxN] do
        for E in EllipticCurves(D,N) do
           T:=A4inftyTorsion(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;

//We build an associative array whose key are possible torsion structure associated
//the corresponding j-maps from table 4. We load the array with the maps that were computed
//in the text and then compute the remaining fiber products.
F<t> := FunctionField(Rationals());
jMaps:=AssociativeArray();
jMaps[[ 1, 1 ]] := <t>;
jMaps[[ 1, 3 ]] := <(t+27)*(t+3)^3/t>;
jMaps[[ 1, 5 ]] := <(t^4-12*t^3+14*t^2+12*t+1)^3/(t^5*(t^2-11*t-1))>;
jMaps[[ 1, 7 ]] := <(t^2+13*t+49)*(t^2+5*t+1)^3/t>;
jMaps[[ 1, 9 ]] := <t^3*(t^3-24)^3/(t^3-27)>;
jMaps[[ 1, 13 ]] := <(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))>;
jMaps[[ 1, 15 ]] := <F!-121945/32, F!46969655/32768>;
jMaps[[ 1, 21 ]] := <F!-140625/8, F!3375/2,F!-1159088625/2097152,F!-189613868625/128>;
jMaps[[ 2, 2 ]] :=  <t^3/(t + 16)>;
jMaps[[ 2, 4 ]] := <(t^6 - 144*t^4 + 6912*t^2 - 110592)/(t^2 - 64)>;
jMaps[[ 2, 6 ]] := <(t+6)^3*(t^3+18*t^2+84*t+24)^3/(t*(t+8)^3*(t+9)^2)>;
jMaps[[ 2, 8 ]] := <(t^12 - 48*t^10 + 816*t^8 - 5632*t^6 + 13056*t^4 - 12288*t^2 + 4096)/(t^4 - 16*t^2)>;
jMaps[[ 2, 10 ]] := <(502270291349/1889568*t^36 - 37469654393/52488*t^35 - 285838091/8748*t^34 + 927749990891/629856*t^33 - 101342580937/209952*t^32 - 50388095567/34992*t^31+ 303116439563/629856*t^30 + 182145354223/209952*t^29 - 294035213/1944*t^28 - 1149479189243/3779136*t^27 - 829060717/139968*t^26 + 806914859/23328*t^25 - 77320765373/16796160*t^24 + 3172839131/699840*t^23 + 1283427409/116640*t^22 + 43829467697/8398080*t^21 + 4181167033/2799360*t^20 + 104705113/155520*t^19 + 3021682909/10077696*t^18 + 17821621/279936*t^17 + 5223229/1866240*t^16 + 281647/16796160*t^15 + 421459/5598720*t^14 - 63853/233280*t^13 - 177396577/1343692800*t^12 - 1258091/111974400*t^11 + 34163/6220800*t^10 + 3236869/2418647040*t^9 - 16553/1343692800*t^8 - 7979/223948800*t^7 - 5507/2015539200*t^6 + 497/1343692800*t^5 + 23/447897600*t^4 - 11/8062156800*t^3 - 1/2687385600*t^2 + 1/967458816000)/(t^36 + 29/4*t^35 + 757/40*t^34 + 1973/160*t^33 - 250331/6400*t^32 - 450947/5120*t^31 - 2428703/64000*t^30 + 4320589/51200*t^29 + 24272243/204800*t^28 + 1572889/81920*t^27 - 28313459/409600*t^26 - 26020517/512000*t^25 + 19052309/3276800*t^24 + 8881181/409600*t^23 + 141671/20480*t^22 - 4873681/1638400*t^21 - 19794679/8192000*t^20 - 194327/819200*t^19 + 856601/3276800*t^18 + 78769/819200*t^17 + 7049/3276800*t^16 - 11957/2048000*t^15 - 2423/1638400*t^14 - 217/1638400*t^13 + 1/409600*t^12 + 1/819200*t^11 + 1/16384000*t^10)>;
jMaps[[ 2, 12 ]] := <(t^24 - 36*t^22 + 522*t^20 - 3924*t^18 + 16551*t^16 - 40392*t^14 + 58860*t^12 - 57672*t^10 + 45279*t^8 - 24084*t^6 + 11178*t^4 - 2916*t^2 + 729)/(t^12 - 12*t^10 + 30*t^8 - 28*t^6 + 9*t^4)>;
jMaps[[ 2, 14 ]] := <F!-3375,F!16581375>;
jMaps[[ 2, 16 ]] := <(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)>;
jMaps[[ 2, 18 ]] := < (t^3-2)^3*(t^9-6*t^6-12*t^3-8)^3/(t^9*(t^3-8)*(t^3+1)^2) >;
jMaps[[ 3, 3 ]] := < 27*(t+1)^3*(t+3)^3*(t^2+3)^3/(t^3*(t^2+3*t+3)^3) >;
jMaps[[ 3, 9 ]] := < F!0 >;
jMaps[[ 4, 4 ]] := <t^2 + 1728, (-64*t^12 + 672*t^10 - 2364*t^8 + 2828*t^6 - 591/4*t^4 + 21/8*t^2 - 1/64)/(t^10 + t^8 + 3/8*t^6 + 1/16*t^4 + 1/256*t^2) >;
jMaps[[ 4, 8 ]] := <(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)>;
jMaps[[ 4, 12 ]] := <(729*t^8 + 756*t^6 + 270*t^4 + 36*t^2 + 1)/t^6>;
jMaps[[ 4, 16 ]] := <(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)>;
jMaps[[ 4, 28 ]] := <(t^4+13*t^2+49)*(t^4+5*t^2+1)^3/t^2>;
jMaps[[ 6, 6 ]]:=<(t^3-57*t^2+3*t-1)^3*(53*t^3+3*t^2-3*t+1)^3*(8587*t^6-8214*t^5+2283*t^4+304*t^3-39*t^2-6*t+1)^3/(729*(t-1)^3*(4*t-1)^6*t^6*(5*t+1)^3*(43*t^2-8*t+1)^3*(7*t^2+t+1)^6)>;
jMaps[[ 8, 8 ]] := <(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)>;



//This function uses the jMaps arrary to compute the E(Q(A4))_tor efficiently.
function A4inftyTorsionFast(E)
    if (Type(E) eq RngIntElt or Type(E) eq FldRatElt) then
        if E eq 1728 or E eq 0 then print "For j-invariant 0 or 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 0 or jE eq 1728 then return A4inftyTorsion(E); end if;
    Tmax:=[1,1];
    for T in Keys(jMaps) do
        for f in jMaps[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(jMaps);
    return Tmax;
end function;

//Does the same thing as A4inftyTorsionExamples but uses A4inftyTorsionFast
function A4inftyTorsionFastExamples(maxN)
    examples:=AssociativeArray();
    D:=CremonaDatabase();
    for N in [11..maxN] do
        for E in EllipticCurves(D,N) do
           T:=A4inftyTorsionFast(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;

//Verifies Table 1 using A4inftyTorsionFase. These will be varified again in ParamCheck2 as long as maxN is set to over 1922.
procedure VerifyTable1()
    print "This will take several minutes...";
    A:=A4inftyTorsionFastExamples(1950);
//    E:=EllipticCurve("1922c1");
//    A[A4inftyTorsion(E)]  := [CremonaReference(E)];
    assert #A eq 26;
    assert A[[ 1, 1 ]][1] eq "11a2";
    assert A[[ 1, 3 ]][1] eq "44a1";
    assert A[[ 1, 5 ]][1] eq "11a1";
    assert A[[ 1, 7 ]][1] eq "26b1";
    assert A[[ 1, 9 ]][1] eq "19a2";
    assert A[[ 1, 13 ]][1] eq "147b1";
    assert A[[ 1, 15 ]][1] eq "50a3";
    assert A[[ 1, 21 ]][1] eq "162b1";
    assert A[[ 2, 2 ]][1] eq "46a1";
    assert A[[ 2, 4 ]][1] eq "17a3";
    assert A[[ 2, 6 ]][1] eq "20a1";
    assert A[[ 2, 8 ]][1] eq "15a5";
    assert A[[ 2, 10 ]][1] eq "66c1";
    assert A[[ 2, 12 ]][1] eq "30a1";
    assert A[[ 2, 14 ]][1] eq "49a1";
    assert A[[ 2, 16 ]][1] eq "210e1";
    assert A[[ 2, 18 ]][1] eq "14a3";
    assert A[[ 3, 3 ]][1] eq "19a1";
    assert A[[ 3, 9 ]][1] eq "27a1";
    assert A[[ 4, 4 ]][1] eq "17a1";
    assert A[[ 4, 8 ]][1] eq "15a2";
    assert A[[ 4, 12 ]][1] eq "30a2";
    assert A[[ 4, 16 ]][1] eq "210e2";
    assert A[[ 4, 28 ]][1] eq "1922c1";
    assert A[[ 6, 6 ]][1] eq "14a1";
    assert A[[ 8, 8 ]][1] eq "15a1";
end procedure;

// Checks the results on CM curves as summarized in Table 3
procedure VerifyTable3()
  assert A4inftyTorsion(EllipticCurve("27a4")) eq [1,9];
  assert A4inftyTorsion(EllipticCurve("32a4")) eq [2,4];
  assert A4inftyTorsion(EllipticCurve("36a2")) eq [2,6];
  assert A4inftyTorsion(EllipticCurve("49a1")) eq [2,14];
  assert A4inftyTorsion(EllipticCurve("121b1")) eq [1,1];
  assert A4inftyTorsion(EllipticCurve("256a1")) eq [2,2];
  assert A4inftyTorsion(EllipticCurve("27a4")) eq [1,9];
  assert A4inftyTorsion(EllipticCurve("361a1")) eq [1,1];
  assert A4inftyTorsion(EllipticCurve("784h2")) eq [2,14];
  assert A4inftyTorsion(EllipticCurve("1849a1")) eq [1,1];
  assert A4inftyTorsion(EllipticCurve("4489a1")) eq [1,1];
  assert A4inftyTorsion(EllipticCurve("26569a1")) eq [1,1];
end procedure;


// Checks the results on CM curves as summarized in Table 4
procedure VerifyTable4()
  assert A4inftyTorsion(EllipticCurve("27a1")) eq [3,9];
  assert A4inftyTorsion(EllipticCurve("36a1")) eq [2,6];
  assert A4inftyTorsion(EllipticCurve("64a4")) eq [4,4];
  assert A4inftyTorsion(EllipticCurve("108a1")) eq [1,3];
  assert A4inftyTorsion(EllipticCurve("256c1")) eq [2,2];
end procedure;

// Checks the last claim in the proof of Prop 4.4
procedure VerifyProp44()
  E:=EllipticCurveFromjInvariant(2268945/128);
  assert pPrimaryA4inftyTorsion(E,7) eq [1,1];
end procedure;

// Checks the last part of the first lemma of the subsubsection The case when p=3
procedure VerifyLemma45()
    G3:=GL(2,Integers(3));
    S0:=[H : H in Subgroups(G3) | IsGenA4(H`subgroup) and FullDeterminantMap(H`subgroup) and ContainsComplexConjugation(H`subgroup)];
    assert #S0 eq 2;
    for H in S0 do
        assert IsConjugateToSubgroup(G3,H`subgroup,GL2SplitCartan(3));
    end for;
end procedure;

// Verifies the first lemma of the subsubsection The case when p=3
procedure VerifyLemma46()
    G9:=GL(2,Integers(9));
    S0:=[H : H in Subgroups(G9) | IsGenA4(H`subgroup) and FullDeterminantMap(H`subgroup) and ContainsComplexConjugation(H`subgroup)];
    for H in S0 do
        assert IsConjugateToSubgroup(G9,H`subgroup,GL2Borel(9));
    end for;
end procedure;



procedure VerifyLemma48()
  S := MaximalImagesWithGenA4Torsion([3,9]);
  assert #S eq 1;
  G := S[1];
  G9 := GL(2,Integers(9));
  G3,mp3 := ChangeRing(G9,Integers(3));
  assert IsConjugateToSubgroup(G9,G,GL2Borel(9));
  assert IsConjugateToSubgroup(G3,mp3(G),GL2Borel(3));
  E:=EllipticCurve("27a1");
  assert A4inftyTorsion(E) eq [3,9];
end procedure;

procedure VerifyRemark49()
  E:=EllipticCurve("54b3");
  assert A4inftyTorsion(E) eq [1,9];
end procedure;


//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 Qab
function TwoPowerLevelStructure(G,n)
    H := ChangeLevel(G,n);
    return GenA4Torsion(H);
end function;


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

//This verifies Corollary 4.11 and requires loading the RZB database
procedure VerifyCorollary411()
n := 8;
printf "Searching the level %o torsion structures of the groups in the RZB database.\n", n;
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[[2,8]];
for i in Tors[[4,8]] do
Append(~Lifts,i);
end for;
for i in Tors[[8,8]] do
Append(~Lifts,i);
end for;
Lifts := Sort(Lifts);
Remove(~Tors,[2,8]);
Remove(~Tors,[4,8]);
Remove(~Tors,[8,8]);

n := 16;
printf "Searching the level %o torsion structures of the groups in the RZB database.\n", n;
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;

Lifts := Tors[[2,16]];
for i in Tors[[4,16]] do
Append(~Lifts,i);
end for;
for i in Tors[[8,16]] do
Append(~Lifts,i);
end for;
assert not [16,16] in Keys(Tors);

Lifts := Sort(Lifts);
Remove(~Tors,[2,16]);
Remove(~Tors,[4,16]);
Remove(~Tors,[8,16]);


assert Lifts eq [ 193, 235, 340, 456, 465, 467, 470, 471, 474, 475, 479, 481, 496, 499, 507, 508, 509, 510, 515, 516, 517, 518, 523, 524, 534, 535, 552, 553, 554, 559, 639,
642, 644, 648, 722, 723, 724, 725 ]; // only 193, 235 have noncuspidal rationals points.

n := 32;
printf "Searching the level %o torsion structures of the groups in the RZB database.\n", n;
for i in [193, 235] 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;
K := Sort([k: k in Keys(Tors)]);
for K in Keys(Tors) do
assert K[2] le 16;
end for;

assert Tors[[ 1, 1 ]] eq [ 1, 3, 4, 5, 7, 20, 22, 55, 56, 57, 177, 178, 179, 180, 439, 440, 441, 718, 719 ];//1 is the maximal representative
assert Tors[[ 2, 2 ]] eq [ 6, 9, 10, 11, 12, 14, 15, 16, 17, 18, 19, 23, 26, 28, 29, 30, 31, 35, 37, 39, 40, 41, 42, 43, 45, 47, 49, 50, 51, 52, 53, 54, 63, 64, 65, 68, 69, 70, 71, 72, 73, 74, 76, 77, 80, 81, 82, 83, 88, 89, 90, 91, 93, 97, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 123, 124, 125, 126, 127, 128, 130, 131, 133, 136, 137, 138, 139, 140, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 160, 161, 162, 163, 164, 165, 166, 167, 169, 170, 171, 172, 173, 174, 175, 176, 191, 196, 216, 218, 220, 224, 232, 237, 238, 239, 250, 252, 253, 257, 258, 261, 262, 263, 264, 265, 266, 267, 275, 276, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 307, 308, 309, 310, 311, 312, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 343, 348, 349, 350, 351, 352, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 370, 372, 373, 374, 375, 377, 380, 382, 383, 384, 385, 387, 388, 389, 390, 391, 394, 396, 397, 399, 400, 401, 402, 403, 414, 415, 416, 417, 418, 419, 420, 421, 422, 426, 443, 461, 480, 484, 520, 525, 526, 533, 538, 539, 542, 545, 546, 547, 548, 549, 550, 551, 555, 556, 557, 558, 563, 564, 565, 566, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 587, 588, 589, 590, 591, 594, 595, 596, 597, 598, 599, 615, 616, 618, 619, 626, 628, 630, 632, 633, 634, 635, 637, 640, 641, 643, 645, 646, 647, 649, 650, 651, 652, 653, 654, 655, 657, 658, 672, 673, 674, 675, 676, 677, 678, 679, 680, 683, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717 ];//6 is the maximal representative
assert Tors[[ 2, 4 ]] eq [ 13, 32, 33, 34, 44, 48, 79, 115, 116, 157, 158, 226, 241, 242, 333, 334, 368, 369, 407, 408, 434, 435, 436, 437, 623, 631, 684, 685 ];//13 is the maximal representative
assert Tors[[ 2, 8 ]] 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 the maximal representative
assert Tors[[ 2, 16 ]] eq [ 235 ];//235 is the maximal representative the rest were removed before the the n = 32 computation becuase they didn't have rational points
assert Tors[[ 4, 4 ]] eq [ 2, 8, 21, 24, 27, 38, 46, 59, 60, 61, 66, 67, 94, 95,132, 134, 135, 141, 209, 210, 214, 221, 222, 231, 254, 256, 260, 277, 376, 378, 379, 381, 392, 393, 395, 398, 411, 412, 413, 425, 537, 543, 600, 601 ]; //2 is the maximal representative
assert Tors[[ 4, 8 ]] 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 the maximal representatives
assert Tors[[ 4, 16 ]] eq [ 193 ]; //193 is the maximal representative the rest were removed before the the n = 32 computation becuase they didn't have rational points
assert Tors[[ 8, 8 ]] 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 the maximal representative
assert not [ 8, 16 ] in Keys(Tors);
end procedure;


procedure Verify413()
  G2:=GL(2,Integers(2));
  G4:=GL(2,Integers(4));
  G:=sub<G2 | [[1,1,0,1]]>;
  invG:=ChangeLevel(G,4);
  S1:=[H`subgroup:H in Subgroups(invG)];
  S2:=[* *];
  for H in S1 do
    if #ChangeLevel(H,2) eq 2 then
      if Exponent(H) eq 2 then
        if IsNormal(G4,H) then
          S2 cat:=[*H*];
        end if;
      end if;
    end if;
  end for;
  assert #S2 eq 0;
end procedure;


procedure VerifyProp53()
    F<t> := FunctionField(Rationals());
	//Construct the generic elliptic curve with a ratonal 13-isogeny. The j-invariant is taken from Table 3 of [18]
    E := EllipticCurveFromjInvariant((t^2+5*t+13)*(t^4+7*t^3+20*t^2+19*t+1)^3/t);
    D := Discriminant(E);
    R := t*(t^2+6*t+13);
    //Verify that D is a square if and only if R is a square
    assert IsSquare(D/R);
    //Define the curve S^2 = R
    E:=EllipticCurve([0,6,0,13,0]);
		assert CremonaReference(E) eq "52a1";
end procedure;

procedure VerifyLemma54()
jlist:=[-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 jlist do
  E:=MinimalTwist(EllipticCurveFromjInvariant(j));
  //print CremonaReference(E);
  assert A4inftyTorsion(E) eq [1,21];
end for;
end procedure;

procedure VerifyLemma55()
F<h> := FunctionField(Rationals());
j := (h^2+13*h+49)*(h^2+5*h+1)^3/h;
E := EllipticCurveFromjInvariant(j);
D := Discriminant(E);
assert IsSquare(D/h);
end procedure;

procedure VerifyRemark56()
  E:=EllipticCurve("1922e2");
  assert A4inftyTorsion(E) eq [4,28];
end procedure;


procedure VerifyLemma57()
  assert A4inftyTorsion(EllipticCurve("49a1")) eq [2,14];
  assert A4inftyTorsion(EllipticCurve("49a1")) eq [2,14];
end procedure;

procedure VerifyLemma59()
  assert A4inftyTorsion(EllipticCurve("50a1")) eq [1,3];
  assert A4inftyTorsion(EllipticCurve("50a2")) eq [1,3];
  assert A4inftyTorsion(EllipticCurve("50a4")) eq [1,15];
  assert A4inftyTorsion(EllipticCurve("50a4")) eq [1,15];
end procedure;


procedure VerifyLemma510()
    F<t> := FunctionField(Rationals());
	//Construct the generic elliptic curve with a ratonal 5-isogeny. The j-invariant is taken from Table 2 of [18]
    E := EllipticCurveFromjInvariant((t^4-12*t^3+14*t^2+12*t+1)^3/(t^5*(t^2-11*t-1)));
    D := Discriminant(E);
  //  DD<x> := Denominator(D);
  //  ND<x> := Numerator(D);
  //  Factorization(DD);
  //  Factorization(ND);
	  R := (t*(t^2-11*t-1));
    //Verify that D is a square if and only if R is a square
    assert IsSquare(D/R);
    //Define the curve S^2 = R
	  E:=EllipticCurve([0,-11,0,-1,0]);
    assert Rank(E) eq 0;
    assert CremonaReference(E) eq "20a4";
end procedure;

procedure VerifyRemark512()
  E:=EllipticCurve("66c1");
  assert A4inftyTorsion(E) eq [2,10];
end procedure;

procedure VerifyRemark515()
  E:=EllipticCurve("30a2");
  assert A4inftyTorsion(E) eq [4,12];
end procedure;


procedure VerifyLemma517()
    F<t> := FunctionField(Rationals());
	//Construct the generic elliptic curve with a rational 3-isogeny. The j-invariant is taken from Table 2 of [18]
    E := EllipticCurveFromjInvariant(((t+27)*(t+3)^3)/t);
    D := Discriminant(E);
    //Verify that D is a square if and only if t is a square
    assert IsSquare(D/t);
end procedure;




/* Here is the code used to compute the various fiber products of the prime power level curves.
We don't run it here but instead sainity check the results in ParamCheck1 and ParamCheck2.

//Compute the j-map for [2,10]
  F<t>:=FunctionField(Rationals());
  j22:=pPrimMap[[2,2]][1];
  j5:=pPrimMap[[1,5]][1];
  R<x,y>:=PolynomialRing(Rationals(),2);
  g:=Numerator(Evaluate(j22,x)-Evaluate(j5,y));
  assert IsIrreducible(g);
  C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
  assert Genus(C) eq 0;
  Co,phi1:=Conic(C);
  P1:=Curve(ProjectiveSpace(Rationals(),1));
  phi2:=Parametrization(Co,P1);
  phi3:=Expand(phi2*Inverse(phi1));
  phi3polys:=DefiningPolynomials(phi3);
  h:=Evaluate(j22,phi3polys[1]/phi3polys[3]);
  jMaps[[2,10]]:=<Evaluate(h,[t,1])>;
//Compute the j-map for [2,12]
  F<t>:=FunctionField(Rationals());
  j24:=pPrimMap[[2,4]][1];
  j3:=pPrimMap[[1,3]][1];
  R<x,y>:=PolynomialRing(Rationals(),2);
  g:=Numerator(Evaluate(j24,x)-Evaluate(j3,y));
  assert IsIrreducible(g);
  C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
  Genus(C);
  Co,phi1:=Conic(C);
  P1:=Curve(ProjectiveSpace(Rationals(),1));
  phi2:=Parametrization(Co,P1);
  phi3:=Expand(phi2*Inverse(phi1));
  phi3polys:=DefiningPolynomials(phi3);
  h:=Evaluate(j24,phi3polys[1]/phi3polys[3]);
  jMaps[[2,12]]:=<Evaluate(h,[t,1])>;
//Compute the j-map for [4,12]
  F<t>:=FunctionField(Rationals());
  assert #pPrimMap[[4,4]] eq 2;
//There are two fiber products to check.
  j441:=pPrimMap[[4,4]][1];
  j442:=pPrimMap[[4,4]][2];
  j3:=pPrimMap[[1,3]][1];
  R<x,y>:=PolynomialRing(Rationals(),2);
  g1:=Numerator(Evaluate(j441,x)-Evaluate(j3,y));
  g2:=Numerator(Evaluate(j442,x)-Evaluate(j3,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 0;
  assert Genus(C2) eq 1;
//Compute the jMaps of the first fiber product that is genus 0
  Co,phi1:=Conic(C1);
  P1:=Curve(ProjectiveSpace(Rationals(),1));
  phi2:=Parametrization(Co,P1);
  phi3:=Expand(phi2*Inverse(phi1));
  phi3polys:=DefiningPolynomials(phi3);
  h:=Evaluate(j441,phi3polys[1]/phi3polys[3]);
  j1:=Evaluate(h,[t,1]);
//Now we show there are no noncuspidal nonsingular rational points on
//the second fiber product that is genus 1
  P := C2![0,0,1];
  E, pi := EllipticCurve(C2,P);
  assert Rank(E) eq 0;
  T , map := TorsionSubgroup(E);
  for P in T do
    assert [R : R in RationalPoints(map(P) @@ pi)] eq [C2![0,0,1],C2![0,1,0],C2![1,0,0]];
  end for;
  jMaps[[4,12]] := <j1>;
//Compute the j-map for [6,6]
  F<t>:=FunctionField(Rationals());
  j22:=pPrimMap[[2,2]][1];
  j33:=pPrimMap[[3,3]][1];
  R<x,y>:=PolynomialRing(Rationals(),2);
  g:=Numerator(Evaluate(j22,x)-Evaluate(j33,y));
  assert IsIrreducible(g);
  C<X,Y,Z>:=ProjectiveClosure(Curve(AffineSpace(R),g));
  assert Genus(C) eq 0;
  Co,phi1:=Conic(C);
  P1:=Curve(ProjectiveSpace(Rationals(),1));
  phi2:=Parametrization(Co,P1);
  phi3:=Expand(phi2*Inverse(phi1));
  phi3polys:=DefiningPolynomials(phi3);
  h:=Evaluate(j22,phi3polys[1]/phi3polys[3]);
  jMaps[[6,6]]:=<Evaluate(h,[t,1])>;

*/



//Here we will not sanity check the parameterization for curves with torsion structure [1,13]
//becuase it is not computationally feasible.
function ParamCheck1(n)
    printf "Sanity checking each parameterization at %o random j-invariants, this will take a while (possibly hours), be patient...\n", n;
    K := Sort([key: key in Keys(jMaps)]);
    for T in K do
        if T ne [1,13] then
        print "Checking parameterization for T =", T;
        for f in jMaps[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:=A4inftyTorsion(E);
                    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 if;
    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)
    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:=A4inftyTorsion(E);
            printf "%o %o\n", CremonaReference(E), T;
            assert A4inftyTorsionFast(E) eq T;
            if jInvariant(E) eq 1728 then assert T in {[2,2],[4,4]}; end if;
            if jInvariant(E) eq 0 then assert T in {[3,9],[2,6], [1,3]}; end if;
        end for;
    end for;
    printf "Successfully checked E(Q(A4^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 3 ***";
    VerifyTable3();
    print "*** Table 3 Verified ***";
    print "*** Verifying Table 4 ***";
    VerifyTable4();
    print "*** Table 4 Verified ***";
    print "*** Verifying Proposition 4.4 ***";
    VerifyProp44();
    print "*** Proposition 4.4 Verified ***";
    print "*** Verifying Lemma 4.5 ***";
    VerifyLemma45();
    print "*** Lemma 4.5 Verified ***";
    print "*** Verifying Lemma 4.6 ***";
    VerifyLemma46();
    print "*** Lemma 4.6 Verified ***";
    print "*** Verifying Lemma 4.8 ***";
    VerifyLemma48();
    print "*** Lemma 4.8 Verified ***";
    print "*** Verifying Remark 4.9 ***";
    VerifyRemark49();
    print "*** Remark 4.9 Verified ***";
    print "*** Verifying Corollary 4.11 ***";
    VerifyCorollary411();
    print "*** Corollary 4.11 Verified ***";
    print "*** Verifying Lemma 4.13 ***";
    Verify413();
    print "*** Lemma 4.13 Verified ***";
    print "*** Verifying Proposition 5.3 ***";
    VerifyProp53();
    print "*** Proposition 5.3 Verified ***";
    print "*** Verifying Lemma 5.4 ***";
    VerifyLemma54();
    print "*** Lemma 5.4 Verified ***";
    print "*** Verifying Lemma 5.5 ***";
    VerifyLemma55();
    print "*** Lemma 5.5 Verified ***";
    print "*** Verifying Remark 5.6 ***";
    VerifyRemark56();
    print "*** Remark 5.6 Verified ***";
    print "*** Verifying Lemma 5.7 ***";
    VerifyLemma57();
    print "*** Lemma 5.7 Verified ***";
    print "*** Verifying Lemma 5.9 ***";
    VerifyLemma59();
    print "*** Lemma 5.9 Verified ***";
    print "*** Verifying Lemma 5.10 ***";
    VerifyLemma510();
    print "*** Lemma 5.10 Verified ***";
    print "*** Verifying Remark 5.12 ***";
    VerifyRemark512();
    print "*** Remark 5.12 Verified ***";
    print "*** Verifying Remark 5.15 ***";
    VerifyRemark515();
    print "*** Remark 5.15 Verified ***";
    print "*** Verifying Lemma 5.17 ***";
    VerifyLemma517();
    print "*** Lemma 5.17 Verified ***";
    print "*** Sanity checking parametrization in Table 4 ***";
    ParamCheck1(3);
    ParamCheck2(2000);
    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 A4inftyTorsionFast(E) is a faster version that uses the parametrizations of Table 4.";
