/////////////////////////////////////////////////////////////////////
// "On the torsion of rational elliptic curves over sextic fields" //
//        Harris B. Daniels and Enrique González-Jiménez           //
/////////////////////////////////////////////////////////////////////

// 2/8/2018 - Magma 2.23

// Magma script related to Proposition 8 (b) 



//load "~/2primary_Ss.txt";
load "2primary_Ss.txt";

//This takes an ordered pair [i,j] and translates it into its corresponding torsion structure. 
//i.e. for rzb in RZB rzb[3][i,j] is the degree of the definition of the (2^k,2^l) torsion,
//this function takes the pair [i,k] and turns it into [2^k,2^l].
function PairToTor(pair)
    if pair[1] eq 1 then return [1,2^pair[2]]; end if;
    return [2^(pair[1]-1),2^(pair[1]+pair[2] -2)];
end function;

//This function takes an item from RZB and returns the corresponding 2-torsion structure over Q. 
function Qtors2(rzb)
    S := Sort([[i,j] : i,j in [1..5]| i+j le 6 and rzb[3][i,j] eq 1]);
    if #S eq 0 then return [1,1]; end if;
    pair := S[#S];
    return PairToTor(pair);
end function;

//Given two torsion groups as ordered pairs, pair1 and pair2, returns true if pair1[1] 
//is divisible by pair2[1] and pair1[2] is divisible by pair2[2]. I.e. if 
// AbelianGroup(pair1) is isomorphic to a subgroup of AbelianGroup(pair2).
function IsSubTor(pair1,pair2)
    return IsDivisibleBy(pair2[1],pair1[1]) and IsDivisibleBy(pair2[2],pair1[2]);
end function;


//Given a list of torsion groups as pairs returns the maximal elements relative to the
//IsSubTor ordering. 
function MaximalTors(list);
    if list eq [] then return [[1,1]]; end if;
    list := Reverse(Sort(list));
    MaxElts := [list[1]];
    for pair in list do
        if {@ false @} eq {@ IsSubTor(pair,L) : L in MaxElts @} then 
            Append(~MaxElts,pair); 
        end if;
    end for;
    return MaxElts;
end function;

//Given an element of RZB returned the an array of torsion structures that occur curves with the 2-adic Galois
//representation base-extended to number fields of degree dividing d. The degree's are the keys of the array.
function Ktors2(rzb, d)
    array := AssociativeArray();
    for fac in Divisors(d) do
    array[fac] := MaximalTors([PairToTor([i,j]) : i,j in [1..5]| i+j le 6 and Integers(rzb[3][i,j])!fac eq 0]);
    end for;
    return array;
end function;

//Takes an array and creates a list of all the items in the array.
function ArrayToList(A)
    L :=[];
    for k in Keys(A) do
    for item in A[k] do
        if not item in L then
            Append(~L,item);
        end if;
    end for;
    end for;
    return L;
end function;

//Searches each element of RZB to what torsion structures appear over degree 6 extensions of Q.
//The results are organized into an associative array whose keys are the torsion over Q.
TorsionArray := AssociativeArray();
for rzb in RZB do
    QT := Qtors2(rzb);
    assert QT in [[1,1],[1,2],[1,4],[1,8],[2,2],[2,4],[2,8] ];
    A := Ktors2(rzb,6);
    L := ArrayToList(A);
    if not QT in Keys(QTorsionArray) then
        QTorsionArray[QT] := L;
    end if;
    if QT in Keys(QTorsionArray) then
    for item in L do
    if not item in QTorsionArray[QT] then
    Append(~QTorsionArray[QT],item);
    end if;
    end for;
    end if;
end for;

assert TorsionArray[[ 1, 1 ]] eq [[ 1 , 1 ],[ 1, 2 ], [ 2, 2 ], [ 1, 4 ], [ 4, 4 ]];
assert TorsionArray[[ 1, 2 ]] eq [[ 1 , 2 ],[ 2, 2 ], [ 1, 4 ], [ 1, 8 ], [ 1, 16 ]];
assert TorsionArray[[ 1, 4 ]] eq [[ 1 , 4 ],[ 2, 4 ], [ 1, 8 ], [ 4, 4 ], [ 2, 8 ]];
assert TorsionArray[[ 1, 8 ]] eq [[ 1 , 8 ],[ 2, 8 ], [ 1, 16 ]];
assert TorsionArray[[ 2, 2 ]] eq [[ 2 , 2 ],[ 2, 4 ], [ 2, 8 ]];
assert TorsionArray[[ 2, 4 ]] eq [[ 2 , 4 ],[ 4, 4 ], [ 2, 8 ]];
assert TorsionArray[[ 2, 8 ]] eq [[ 2 , 8 ]];

//Got <----

// 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 a subgroup of GL2(Z/mZ) lift or reduce it a subgroup of GL2(Z/nZ)
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;


// 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 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 list of groups of GL(2,Integers(N)) return a list of the maximal elements.
function MaximalGroups(S,N)
G := GL(2,Integers(N));
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;
    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;



function DegreeDTorStructures2Primary(QT,KT,degree)
    DATA := <>;
    for rzb in RZB do
        if not QT in [[1,1],[1,2],[1,4],[1,8],[2,2],[2,4],[2,8] ] then return []; end if;
        if Qtors2(rzb) eq QT then
        A := Ktors2(rzb,degree);
        L := ArrayToList(A);
        if KT in L then
        Append(~DATA,rzb);
        end if;
        end if;
    end for;
    if #DATA eq 0 then return []; end if;
    N := GCD([data[2] : data in DATA]);
    Groups := [sub<GL(2,Integers(d[2]))|[a : a in d[4]]> : d in DATA ];
    Labels := [d[1] : d in DATA ];
    Lifts := [ChangeLevel(G,N) : G in Groups];
    Max := MaximalGroups(Lifts,N);
    MaxCurves := [rzb[1] : rzb in RZB | sub<GL(2,Integers(rzb[2]))|[a : a in rzb[4]]> in Max];
    assert #Max eq #MaxCurves;
    return MaxCurves;
end function;


//This function takes in a degree and returns an associative array whose keys are 
//torsion configurations [G,H] and the corresponding object is a sequence of labels
//modular curve parameterizing elliptic curves with rational torsion G and have torsion 
//structure H over a degree deg extension of Q.
function TwoTorsConfiguration(deg)
    QTorsionArray := AssociativeArray();
    for rzb in RZB do
        QT := Qtors2(rzb);
        assert QT in [[1,1],[1,2],[1,4],[1,8],[2,2],[2,4],[2,8] ];
        A := Ktors2(rzb,deg);
        L := ArrayToList(A);
        if not QT in Keys(QTorsionArray) then
            QTorsionArray[QT] := L;
        end if;
        if QT in Keys(QTorsionArray) then
        for item in L do
        if not item in QTorsionArray[QT] then
            Append(~QTorsionArray[QT],item);
        end if;
        end for;
        end if;
    end for;
    ModularCurvesArray := AssociativeArray();
    for QT in Keys(QTorsionArray) do
    for KT in QTorsionArray[QT] do
        ModularCurvesArray[[QT,KT]] := DegreeDTorStructures2Primary(QT,KT,6);
    end for;
    end for;
    return ModularCurvesArray;
end function;

A := TwoTorsConfiguration(6);

assert A[[[ 1, 1 ],[ 1, 1 ]]] eq  [ "X1" ];
assert A[[[ 1, 1 ],[ 1, 2 ]]] eq  [ "X1" ];
assert A[[[ 1, 1 ],[ 1, 4 ]]] eq  [ "X20" ];
assert A[[[ 1, 1 ],[ 2, 2 ]]] eq  [ "X1" ];
assert A[[[ 1, 1 ],[ 4, 4 ]]] eq  [ "X20b" ];
assert A[[[ 1, 2 ],[ 1, 2 ]]] eq  [ "X6" ];
assert A[[[ 1, 2 ],[ 1, 4 ]]] eq  [ "X13" ];
assert A[[[ 1, 2 ],[ 1, 8 ]]] eq  [ "X102", "X36a" ];
assert A[[[ 1, 2 ],[ 1, 16 ]]] eq  [ "X235m" ];
assert A[[[ 1, 2 ],[ 2, 2 ]]] eq  [ "X6" ];
assert A[[[ 1, 4 ],[ 1, 4 ]]] eq  [ "X13h" ];
assert A[[[ 1, 4 ],[ 1, 8 ]]] eq  [ "X36n" ];
assert A[[[ 1, 4 ],[ 2, 4 ]]] eq  [ "X13h" ];
assert A[[[ 1, 4 ],[ 2, 8 ]]] eq  [ "X102k" ];
assert A[[[ 1, 4 ],[ 4, 4 ]]] eq  [ "X60d" ];
assert A[[[ 1, 8 ],[ 1, 8 ]]] eq  [ "X102p" ];
assert A[[[ 1, 8 ],[ 1, 16 ]]] eq  [ "X235l" ];
assert A[[[ 1, 8 ],[ 2, 8 ]]] eq  [ "X102p" ];
assert A[[[ 2, 2 ],[ 2, 2 ]]] eq  [ "X8" ];
assert A[[[ 2, 2 ],[ 2, 4 ]]] eq  [ "X25", "X8d" ];
assert A[[[ 2, 2 ],[ 2, 8 ]]] eq  [ "X193", "X96q", "X98o" ];
assert A[[[ 2, 4 ],[ 2, 4 ]]] eq  [ "X25n" ];
assert A[[[ 2, 4 ],[ 2, 8 ]]] eq  [ "X96t", "X98e" ];
assert A[[[ 2, 4 ],[ 4, 4 ]]] eq  [ "X58i" ];
assert A[[[ 2, 8 ],[ 2, 8 ]]] eq  [ "X193n" ];