///////////////////////////////////////////////////////////////////////////
//
// 		"An algorithm for determining torsion growth of elliptic curves"
// 		Enrique González-Jiménez & Filip Najman
//
// 		27/3/2019 - Magma 2.24
//      file: programNF_optimized2.m
//
///////////////////////////////////////////////////////////////////////////
booprint:=false;
booprint2:=false;
booprint3:=false;
booprintGal:=true;
booOptimized:=false;
dOR:=16; // maximun degree to use OptimizedRepresentation
printf "-----------------------------------------------------\n Printing level (1,2,3,Gal)=(%o,%o,%o,%o)\n",booprint,booprint2,booprint3,booprintGal;
if booOptimized then 
	printf " Using Optimized Representation for degree <=%o? %o\n-----------------------------------------------------\n",dOR,booOptimized; 
  else 
    printf " Using Optimized Representation? %o\n-----------------------------------------------------\n",booOptimized; 
end if;
///////////////////////////////////////////////////////////////////////////

load "step1NF.m";
load "step2NF_optimized.m";
load "step3NF.m";

///////////////////////////////////////////////////////////////////////////

function OrderGalrep(E, ell)
local a,b;
if not HasComplexMultiplication(E) then
	galrep:=GaloisLabel(E,ell);
	if Substring(galrep,1,3) eq "GL2" then
		return ell*(ell^2-1)*(ell-1);
	else
		a,b:=ExcepcionalGalrepDegrees(galrep);
		return b;
	end if;
else
	if ell eq 2 then a,b:=DegCM_2D(E);
	else a,b:=DegreesCM(E, ell);
	end if;
return b;
end if;
end function;

//two auxiliary functions
function Degrees(A)
return [Degree(b[1]):b in A];
end function;

function aux_Groups(A)
return <[bb : bb in b[2]]:b in A>;
end function;


function TorsionGrowth(E,d)
if booprint then printf "%o\n", "//// TORSION GROWTH ////"; T0 := Cputime(); end if;
E:=MinimalModel(E);
lst:=<>;
primes:=R(E,d);
if booprint then printf "R(E,%o)=%o\n",d,primes; end if;
for i:=1 to #primes do
	ell:=primes[i];
	if booprint then printf "*** about to compute step2 for ell=%o ***\n",ell; t0 := Cputime(); end if;
	A:=step2(E,ell,d, OrderGalrep(E, ell));
	if booprint then printf "*** step2 for ell=%o done in %o seconds ***\n",ell,Cputime(t0); end if;
	lst:=lst cat A;
end for;



T:=TorsionSubgroup(E);
fact:=Factorization(#T);
for i:=1 to #fact do
	ell:=fact[i,1];
	if ell notin primes then
	g:=pPowerTorsion(E,ell);
	if Exponent(g) eq Order (g) then
		lst:= lst cat <<RationalsAsNumberField(), <Exponent(g)>, {<ell, RationalsAsNumberField()>}>>;
	else 
	lst:= lst cat <<RationalsAsNumberField(), <Exponent(g), Order(g) div Exponent(g)>, {<ell, RationalsAsNumberField()>}>>;
	end if;
	end if;
end for;

if booprint then printf "****** %o...", "about to perform step3"; t0 := Cputime(); end if;
B:=step3(lst,d);
if booprint then printf "****** step 3 done in %o seconds\n",Cputime(t0); end if;

A:=<<NumberField(Polynomial(a[1])),<aa : aa in a[2]>> : a in B>;
B:=<>; 

//Degrees(A);
//aux_Groups(A);
for i:=1 to #A do
	check:=true;
	for j:= 1 to #B do
	if booprint then printf "[i,j]=[%o,%o]\n", i,j; end if;
	if booprint then printf "%o","about to test IsIsomorphicNF..."; t0 := Cputime(); end if;	
	Ai:=A[i,1]; Bj:=B[j,1]; 
	if booprint then printf "\n%o\n%o\n", Ai,Bj; end if;
	booij:=IsIsomorphicNF(Ai,Bj);
	if booprint then printf "done in %o seconds\n",Cputime(t0); end if;	
	if booij then
		//i,j,A[i,2], B[j,2];
		check:=false;
		if A[i,2,1] ge B[j,2,1] then
			if (#B[j,2] eq 1) then
				C:=<>;
				for numer:=1 to #B do
					if numer ne j then
						C:=C cat <B[numer]>;
					else 
						assert numer eq j;
						C:=C cat <A[i]>;
					end if;
				end for;
				B:=C;
			else
				if #A[i,2] eq 2 then
					if A[i,2,2] ge B[j,2,2] then
						B[j,1]:=A[i,1]; B[j,2]:=A[i,2];
					end if;
				end if;
			end if;
		else
			check:=false;		
		end if;
	end if;
	end for;
	//print check;
	if check then 
		B:= B cat <A[i]>;
	end if;
end for;
if booprint then printf "//// DONE IN %o seconds ////\n", Cputime(T0); end if;
if booOptimized then B:=<<SwitchToOptimized(K[1]),K[2]> : K in B>; end if;
return [[[aa : aa in a[2]],Coefficients(DefiningPolynomial(a[1]))] : a in B];
end function;