///////////////////////////////////////////////////////////////////////////
//
// 		"An algorithm for determining torsion growth of elliptic curves"
// 		Enrique González-Jiménez & Filip Najman
//
// 		6/3/2019 - Magma 2.24
//      file: step3NF.m
//
///////////////////////////////////////////////////////////////////////////
//booprint3:=false;
booCompoNF:=true;
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// We suppost that we receive as input a list A consisting of pairs (K,T), where K is a number field and T is an l-torison group, 
// represented as either one ell-power (if the group is cyclic).
function step3(A,d)
// local  t;
B:=<>;
for i:= 1 to #A do
	ell:=Factorization(A[i,2,1])[1,1];
	temppair:=A[i][2];
	if #temppair eq 2 then
		temppgroup:=DirectProduct(CyclicGroup(temppair[1]), CyclicGroup(temppair[2]));
	else 
		assert #temppair eq 1;
		temppgroup:=CyclicGroup(temppair[1]);
	end if;
	temp:=<<A[i][1],temppgroup, { <ell, A[i][1]> } >>;
	B:=B cat temp;
end for;


k:=2;
repeat
	if booprint3 then printf "k=%o\n", k; end if;
	new:=false;
	nmr:=#B;
	for i:=1 to nmr do
		for j:=i+1 to nmr do
			if booprint3 then printf "[i,j]=[%o,%o]\n", i,j; end if;
			Sp:={a[1]:a in B[i][3]} join  {a[1]:a in B[j][3]};
			if (#Sp eq k) and (#(B[i][3] join B[j][3]) eq k) then
				//i,j, Degree(B[i,1]), Degree(B[j,1]), #B[i,2], #B[j,2];
				if Degree(B[i][1]) le Degree(B[j,1]) then 
					K:=B[i][1]; p:=(PolynomialRing(K))!DefiningPolynomial(B[j][1]);
				else 
					K:=B[j][1]; p:=(PolynomialRing(K))!DefiningPolynomial(B[i][1]);
				end if;
				d1:=d div Degree(K);
				test:=false;
				if booprint3 then printf "%o","about to test IsSubfield..."; t0 := Cputime(); end if;	
				booij:=IsSubfieldNF(B[j,1], B[i,1]);
				if booprint3 then printf "done in %o seconds\n",Cputime(t0); end if;				
				if Degree(B[i,1]) eq d and (not booij) then 
					test:=true; 
				else 
					if booprint3 then printf "%o","about to test IsSubfield..."; t0 := Cputime(); end if;	
					booij:=IsSubfieldNF(B[i,1], B[j,1]);
					if booprint3 then printf "done in %o seconds\n",Cputime(t0); end if;	
					if Degree(B[j,1]) eq d and (not booij) 
					then test:=true; 
					else 
						//test:=poltest(p,&*Sp, d1);
						if booprint3 then printf "%o","about to factor polynomial..."; t0 := Cputime(); end if;
						fac:=Factorization(p);
						if booprint3 then printf "done in %o seconds\n",Cputime(t0); end if;
						d1:= d div Degree(K);
						if d1 mod Degree(fac[1][1]) eq 0 then 
							test:= false;
						else
							test:=true;
						end if;
					end if;
					//print "done";
				end if;
				if test eq false then
					if booprint3 then printf "%o","about to compute Compositum..."; t0 := Cputime(); end if;
					if booCompoNF then 
						if booprint3 then printf "%o","using NF version..."; end if;
						D:=CompositumNF(B[i][1], B[j][1]); 
					  else
					    if booprint3 then printf "%o","using classical version..."; end if;
					    D:=CompositeFields(B[i][1], B[j][1])[1];
					end if;
					if booprint3 then printf "done in %o seconds\n",Cputime(t0); end if;
					//#CompositeFields(B[i][1], B[j][1]),B[i][1], B[j][1], D;
					dp:=Degree(D);
					if (d mod dp) eq 0 then
						new:=true;
						//Sp,k, #B[i][2],#B[j][2], B[j][3] diff B[i][3];
						T:=DirectProduct(B[i][2], DirectProduct([SylowSubgroup(B[j][2],a[1]):a in (B[j][3] diff B[i][3])]));
						//#T, B[i][3] join B[j][3];
						B:=B cat <<D, T, B[i][3] join B[j][3]>>;
						//#T;
					end if;
				end if;
			end if;
		end for;
	end for;
	k:=k+1;
until new eq false;
C:=<>;
for i:=1 to #B do
tempgroup:=B[i,2];
if IsCyclic(tempgroup) then
		temppair:=<Order(tempgroup)>;
	else 
		exp:=Exponent(tempgroup);
		ord:=Order(tempgroup);
		temppair:=<exp, ord div exp>;
	end if;	
	//Degree(B[i,1]), temppair;
	C:=C cat <<B[i,1], temppair>>;
end for;
return [[Coefficients(DefiningPolynomial(a[1])),[aa : aa in a[2]]] : a in C];
end function;

