/*
	This code is for the paper "PUT THE TITLE HERE" by Harris B. Daniels and Alvaro Lozano-robledo

	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)
*/




/*
This code requires loading the databases of groups from [18].
This data is available at

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

*/

//Load the databases from RZB
load "gl2data";
load "gl2finedata";

print "Run VerifyAll() in order to verify all the computations";


F<t> := FunctionField(Rationals());
//Function that returns the jmap for a modular curve from the RZB database.
function jmap(subnum)
  curmap := F!t;
  done := false;
  curgp := subnum;
  while (done eq false) do
    if (curgp eq 1) then
      done := true;
      j := curmap;
    else
      covernum := newsublist[curgp][7];
      filename := Sprintf("%omap%o.txt",curgp,covernum);
      filename1 := Sprintf("%map%o.txt",curgp,covernum);
      maptocover := eval Read(filename);
      newcurmap := maptocover([ curmap, 1])[1];
      curmap := newcurmap;
      curgp := covernum;
    end if;
  end while;
  return curmap;
end function;

// 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;



// 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;

// 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;


// 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 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;


//Change the level of G a subgroup of GL(2,Z/mZ).
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;

//Takes the transpose of a subgroup of a matrix group.
function MatrixGroupTranspose(H)
	G := GL(2,BaseRing(H));
    return sub<G| [Transpose(A) : A in Generators(H) ]>;
end function;


/* 
Example 1.2
*/
procedure VerifyEx12()
	E := EllipticCurve([405,-9882]);
	assert GroupName(TorsionSubgroup(E)) eq "C3";
	f2 := DivisionPolynomial(E,2);
	K<a> := SplittingField(f2);
	EK := BaseExtend(E,K);
	assert GroupName(TorsionSubgroup(EK)) eq "C6^2";
end procedure;

/* 
Example 1.3
*/
procedure VerifyEx13()
	E := EllipticCurve([13,-34]);
	assert GroupName(TorsionSubgroup(E)) eq "C4";
	f2 := DivisionPolynomial(E,2);
	K<i> := NumberField(Polynomial([1,0,1]));
	EK := BaseExtend(E,K);
	assert GroupName(TorsionSubgroup(EK)) eq "C4^2";
	P := EK![7,-20];
	Q := EK![-3-10*i,30+10*i];
	assert Order(P) eq 4 and Order(Q) eq 4;
	assert EK![-3+10*i,30-10*i] eq P+3*Q;
end procedure;

/*
Verify the family in Theorem 1.6
*/

procedure VerifyFamily16()
	_<x> := PolynomialRing(Rationals());
	K<sqrt> := NumberField(x^2+3);
	F<t> := FunctionField(K);
	A := -3*(-(t^3+1)^2+1)^3*(-(t^3+1)^2+9);
	B := -2*(-(t^3+1)^2+1)^4*((-(t^3+1)^2)^2-18*(-(t^3+1)^2)-27);
	E := EllipticCurve([ A, B ]);
	f2 := DivisionPolynomial(E,2);
	F2 := ext<F|f2>;
	E2 := BaseExtend(E,F2);
	roo2 := Roots(DivisionPolynomial(E2,2));
	assert #roo2 eq 3;
	roo3 := Roots(DivisionPolynomial(E2,3));
	a1 := roo3[1][1];
	a2 := roo3[2][1];
	P<y> := PolynomialRing(F2);
	_,b1 := IsSquare(a1^3+A*a1+B);
	_,b2 := IsSquare(a2^3+A*a2+B);
	P1 := (E2![a1,b1]);
	P2 := (E2![a2,b2]);
	assert Order(P1) eq 3 and Order(P2) eq 3;
	assert P1 ne P2 and P1 ne 2*P2;
end procedure;


/* 
Example 3.7
*/
procedure VerifyEx37()
	E := EllipticCurve([1,-1,1,4,-1]);
	assert GroupName(pPowerTorsion(E,2)) eq "C1";
	f2:=DivisionPolynomial(E,2);
	assert GroupName(GaloisGroup(f2)) eq "S3";
	K<a> := SplittingField(f2);
	assert DefiningPolynomial(K) eq Polynomial([82944,0,42849,0,6624,0,256]);
	EK := BaseExtend(E,K);
	assert GroupName(pPowerTorsion(EK,2)) eq "C4^2";
end procedure;

procedure VerifyProp38()
	F<t> := FunctionField(Rationals());
	G4 := GL(2,Integers(4));
	G2, mp2 := ChangeRing(G4,Integers(2));
	SG := Subgroups(G4);

	// Get all the subgroups H of GL(2,Integers(4)) that mp2(H) \simeq H.
	list := [h`subgroup : h in SG | FullDeterminantMap(h`subgroup) and #(h`subgroup meet Kernel(mp2)) eq 1];

	// All the CM j-invariants
	CMj := [ 0, 54000, -12288000, 1728, 287496, -3375, 16581375, 8000, -32768, -884736,
	-884736000, -147197952000, -262537412640768000 ];

	//Showing all the possible images are contained in the group X_{20b} which we call Im in the code.
	Im := MatrixGroupTranspose(finesublist[26][4]);

	for H in list do
		assert IsConjugateToSubgroup(GL(2,Integers(4)), H, Im);
	end for;


	//Define jmp to be the jmap for X_{20}
    jmp := (-4*t^8 + 32*t^7 + 80*t^6 - 288*t^5 - 504*t^4 + 864*t^3 + 1296*t^2 - 864*t - 1188)/(t^4 + 4*t^3 + 6*t^2 + 4*t + 1);

	//Show that non of the CM j-invariants are in the image of jmp.
	for j in CMj do
		rts := Roots(Numerator(jmp - j));
		assert #rts eq 0;
	end for;

	//Set up
	P<t>:=FunctionField(Rationals());
	A := -27*t^(8) + 648*t^(7) - 4212*t^(6) - 2376*t^(5) + 60102*t^(4) + 79704*t^(3) - 105732*t^(2) - 235224*t - 107811;
	B := 54*t^(12) - 1944*t^(11) + 24300*t^(10) - 97848*t^(9) - 251262*t^(8) + 1722384*t^(7) + 4821768*t^(6) - 8697456*t^(5) - 64323558*t^(4) - 140447736*t^(3) - 157012020*t^(2) - 90561240*t - 21346578;
	E:=EllipticCurve([A,B]);
	assert jInvariant(E) eq jmp;

	//Code to verify Q(E[2^n])=/= Q(E[2^{n+1}]) for j=1728

	n:=32;
	Gn:=GL(2,Integers(n));

	d:=-4;
	L:=[[a,b,-b,a] : a, b in [0..(n-1)] | ((a^2+b^2) mod 2) eq 1];
	L:=Append(L,[-1,0,0,1]);
	G4a:=sub<Gn|L>;

	d:=-4;
	G4b:=sub<Gn|[[-1,0,0,1],[-1,0,0,-1],[3,0,0,3],[1,2,-2,1]]>;

	d:=-4;
	G4b2:=sub<Gn|[[0,-1,-1,0],[-1,0,0,-1],[3,0,0,3],[1,2,-2,1]]>;


	d:=-4;
	G4c:=sub<Gn|[[0,-1,-1,0],[-1,0,0,-1],[3,0,0,3],[2,1,-1,2]]>;

	d:=-4;
	G4c2:=sub<Gn|[[-1,0,0,1],[-1,0,0,-1],[3,0,0,3],[2,1,-1,2]]>;


	d:=-4;
	G4d:=sub<Gn|[[-1,0,0,1],[5,0,0,5],[1,2,-2,1]]>;

	d:=-4;
	G4d2:=sub<Gn|[[0,-1,-1,0],[5,0,0,5],[1,2,-2,1]]>;

	d:=-4;
	G4e:=sub<Gn|[[-1,0,0,1],[5,0,0,5],[-1,-2,2,-1]]>;

	d:=-4;
	G4e2:=sub<Gn|[[0,-1,-1,0],[5,0,0,5],[-1,-2,2,-1]]>;


	d:=-4;
	G4f:=sub<Gn|[[1,0,0,-1],[-3,0,0,-3],[2,-1,1,2]]>;

	d:=-4;
	G4f2:=sub<Gn|[[-1,0,0,1],[-3,0,0,-3],[2,-1,1,2]]>;


	d:=-4;
	G4g:=sub<Gn|[[1,0,0,-1],[-3,0,0,-3],[-2,1,-1,-2]]>;

	d:=-4;
	G4g2:=sub<Gn|[[-1,0,0,1],[-3,0,0,-3],[-2,1,-1,-2]]>;

	images:=[G4a,G4b,G4b2,G4c,G4c2,G4d,G4d2,G4e,G4e2,G4f,G4f2,G4g,G4g2];
	for G in images do
		assert #ChangeRing(G,Integers(4)) div #ChangeRing(G,Integers(2)) ge 2 and #ChangeRing(G,Integers(8)) div #ChangeRing(G,Integers(4)) ge 2 and #ChangeRing(G,Integers(16)) div #ChangeRing(G,Integers(8)) ge 2  and #ChangeRing(G,Integers(32)) div #ChangeRing(G,Integers(16)) ge 2;
	end for;

end procedure;


/* 
Search the finesublist
*/
procedure VerifySection4()

for array in finesublist do 
	level := array[3];
	exponent := Factorization(level)[1][2];
	group := array[4];
	gens := Generators(group);
	if exponent ne 1 then
		for n in [1..exponent-1] do
		G1 := sub<GL(2,Integers(2^n))|gens>;
		G2 := sub<GL(2,Integers(2^(n+1)))|gens>;
		if IsIsomorphic(G1,G2) then
			assert array[5] in { "X20b", "X60d" };
		end if;
		end for;
    end if;
end for;


for i in [1..#newsublist] do
	group := newsublist[i][3];
    level := #BaseRing(group);
	exponent := Factorization(level)[1][2];
	gens := Generators(group);
	if exponent ne 1 then
		for n in [1..exponent-1] do
		G1 := sub<GL(2,Integers(2^n))|gens>;
		G2 := sub<GL(2,Integers(2^(n+1)))|gens>;
		if IsIsomorphic(G1,G2) then
			assert 1 eq 0;
		end if;
		end for;
    end if;
end for;

end procedure;



/*
Theorem 1.7
*/
procedure VerifyTheorem17()
	G6 := GL(2,Integers(6));
    SG := [ h : h in Subgroups(G6) | FullDeterminantMap(h`subgroup) and ContainsComplexConjugation(h`subgroup) ];
    list := [h`subgroup : h in SG | IsIsomorphic(h`subgroup, ChangeLevel(h`subgroup,2)) and IsIsomorphic(h`subgroup, ChangeLevel(h`subgroup,2)) and IsAbelian(h`subgroup)];
    for h in list do
		assert IsConjugate(G6, h, sub<G6| [1,1,0,5]>);
    end for;
end procedure;



/*
Lemma 6.2
*/
procedure VerifyLemma62()

G := GL(2,Integers(12));
SG := [h`subgroup : h in Subgroups(G) | FullDeterminantMap(h`subgroup)];

test := [];
test3 := [];
test4 := [];
for H in SG do
	H4 := ChangeLevel(H,4);
	H3 := ChangeLevel(H,3);
	if not IsAbelian(H) and IsIsomorphic(H3,H4) and IsIsomorphic(H,H3) then
		Append(~test,H);
		Append(~test3,H3);
		Append(~test4,H4);
	end if;
end for;

assert #test eq 2;

for H in test4 do
	assert GroupName(H) eq "D4";
end for;

assert IsConjugate(GL(2,Integers(4)), test4[2], finesublist[14][4]);
assert IsConjugate(GL(2,Integers(4)), test4[1], finesublist[16][4]);

end procedure;


/*
Lemma 6.4
*/

procedure VerifyLemma64()
G9 := GL(2,Integers(9));
SG9 := [h`subgroup : h in Subgroups(G9) | not IsAbelian(h`subgroup) and FullDeterminantMap(h`subgroup) and ContainsComplexConjugation(h`subgroup)];
S := { "C6", "C2*C6", "C3*C6", "C6^2" };
for H in SG9 do
	assert GroupName(H/DerivedGroup(H)) in S;
end for;



H31 := sub<GL(2,Integers(3)) |  [1,0,0,2],[2,1,2,2] >; //N_{ns}(3)
H32 := GL2Borel(3);
H91 := ChangeLevel(H31,9);
H92 := ChangeLevel(H32,9);

S1 := [H : H in SG9 | GroupName(H/DerivedGroup(H)) eq "C2*C6" ];

for K in S1 do
	assert IsConjugateToSubgroup(GL(2,Integers(9)),K, H91) or IsConjugateToSubgroup(GL(2,Integers(9)),K, H92);
end for;


G4 := GL(2,Integers(4));
SG4 := [h`subgroup : h in Subgroups(G4) | not IsAbelian(h`subgroup) and FullDeterminantMap(h`subgroup) and ContainsComplexConjugation(h`subgroup)];
T := { "C2*C4", "C2^3", "C2*C6", "C2^2", "C6", "C2" };

assert {GroupName(h/DerivedGroup(h)) : h in SG4 } eq T;

S2 := [H : H in SG4 | GroupName(H/DerivedGroup(H)) eq "C2*C6" ];

H2 := sub<GL(2,Integers(2)) | [0,1,1,1] >;
H4 := ChangeLevel(H2,4);

for K in S2 do
	assert IsConjugateToSubgroup(GL(2,Integers(4)),K, H4) ;
end for;



Subs91 := [ h`subgroup : h in Subgroups(H91) | FullDeterminantMap(h`subgroup) and ContainsComplexConjugation(h`subgroup) ];
Subs92 := [ h`subgroup : h in Subgroups(H92) | FullDeterminantMap(h`subgroup) and ContainsComplexConjugation(h`subgroup) ];
Subs4 := [ h`subgroup : h in Subgroups(H4) | FullDeterminantMap(h`subgroup) and ContainsComplexConjugation(h`subgroup) ];



for K4 in Subs4 do
if true in {IsIsomorphic(K9,K4) : K9 in Subs91} then
    assert 1 eq 0;
end if;
end for;

for K4 in Subs4 do
if true in {IsIsomorphic(K9,K4) : K9 in Subs92} then
    assert 1 eq 0;
end if;
end for;

end procedure;


procedure VerifyAll()
	print "Verifying all computations... This may take a few minutes...";
    print " ";
	print "Verifying Example 1.2";
    VerifyEx12();
    print "Example 1.2 Verified";
    print " ";
    print " ";
    print "Verifying Example 1.3";
    VerifyEx13();
    print "Example 1.3 Verified";
    print " ";
    print " ";
    print "Verifying the Family in Theorem 1.6";
    VerifyFamily16();
    print "Verified the Family in Theorem 1.6";
    print " ";
    print " ";
    print "Verifying Example 3.7";
    VerifyEx37();
	print "Verified Example 3.7";
    print " ";
    print " ";
    print "Verifying Proposition 3.8";
    VerifyProp38();
    print "Verified Proposition 3.8";
    print " ";
    print " ";
    print "Verifying Section 4 Search";
    VerifySection4();
    print "Verified the Search in Section 4";
    print " ";
    print " ";
    print "Verifying Theorem 1.7";
    VerifyTheorem17();
    print "Verified Theorem 1.7";
    print " ";
    print " ";
    print "Verifying Lemma 6.2";
    VerifyLemma62();
    print "Verified Lemma 6.2";
    print " ";
    print " ";
    print "Verifying Lemma 6.4";
    VerifyLemma64();
    print "Verified Lemma 6.4";
    print " ";
    print " ";
    print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
    print "%%%%% All The Computations Have Been Verified. %%%%%";
    print "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
end procedure;




