///////////////////////////////////////////////////////////////////////////
//
// 		"An algorithm for determining torsion growth of elliptic curves"
// 		Enrique González-Jiménez & Filip Najman
//
// 		27/3/2019 - Magma 2.24
//      last modified: 13/5/2020
//      file: step2NF_optimized.m
//
///////////////////////////////////////////////////////////////////////////

//booprint2:=false;
//booOptimized:=true;
// dOR:=8; // maximun degree to use OptimizedRepresentation
// the function below takes a E/K and a number a in K and returns the field over which the point on E(K) with a coordinate is defined, together with a point
// E and a should be defined over the same field

function OrderC(P)
	zfunc:=0;
	repeat zfunc:=zfunc+1; until IsZero(zfunc*P); 
	return zfunc;
end function;

function SwitchToOptimized(K)
	if Degree(K) le dOR and booOptimized then 
		F,f:=OptimizedRepresentation(K);
		return F,f;
	  else 
		return K,IdentityFieldMorphism(K);
	end if;
end function;



function NFext(E,a)
// local K,a1,a2,a3,a4,a6,y2,booIrre,L,tr,f,f2,L1, P, pts,E1;
assert IsIsomorphicNF(BaseField(E), Parent(a));
K:=BaseField(E);
a1,a2,a3,a4,a6:=Explode(aInvariants(E));
y2:=Polynomial([-(a^3+a2*a^2+a4*a+a6),(a1*a+a3),1]);
booIrre:=IsIrreducible(y2);
if booIrre then
	L:=ext<K|y2>;
	//if booprint then printf "%o","about to compute isomorphism..."; t0 := Cputime(); end if;
	//tr,f:=IsSubfieldNF(K,L);
	//assert tr;
	L1:=AbsoluteField(L);
	//polynom:=AbsoluteMinimalPolynomial(a);
	//assert not IsIrreducible(PolynomialRing(L1)!polynom);
	//atemp:=Roots(PolynomialRing(L1)!polynom)[1,1];
	//tr,f2:=IsIsomorphic(L,L1);
	//if booprint then printf "done in %o seconds\n",Cputime(t0); end if;
	//assert tr;
	//atemp:=f2(L!a);
	atemp:=L1!a;
	E1:=ChangeRing(E,L1);
	pts:=Points(E1,atemp);
	assert #pts ge 1;
	P:=pts[1];
	return L1,P;
else
	L1:=AbsoluteField(K);
	pts:=Points(ChangeRing(E,L1),L1!a);
	assert #pts ge 1;
	P:=pts[1];
	return L1,P;
end if;
end function;

// the function below replaces a triple with a different triple, it is neccessary this way as magma won't allow in an easier way

function replace(A,i,t)
// local numer, C;
C:=<>;
	for numer:=1 to #A do
		if numer ne i then
			C:=C cat <A[numer]>;
		else 
			assert numer eq i;
			C:=C cat <t>;
		end if;
	end for;
return C;				
end function;


function DivisionFieldEasy(E,n)
// E an elliptic curve over a number field K
// n an odd integer
// we assume that there is a point of order n on E defined over K, so now it's enough to add all \zeta_n and 
// all the x-coordinates of the n-division polynomaial (by Bandini-Paladino, JNT, Theorem 1.1.)
//
// local K,K1,i,lst,fac,T,m,E1,facx0,P0,S,z0,P1;
p:=PrimitiveDivisionPolynomial(E,n);
fac:=Factorization(p);
//pos:=false;
i:=0;
K:=BaseField(E);
K1,lst:=SplittingField(p);
_<x>:=PolynomialRing(K1);
fac2:=Factorization(x^n-1);
K2:=ext<K1|fac2[#fac2][1]>;
K:=AbsoluteField(K2);
T,m:=pPowerTorsion(E,n);
Ptemp:=m(T.1);
mult:= OrderC(Ptemp) div n;
P0:=mult*Ptemp;
assert OrderC(P0) eq n;
if OrderC(P0) ne n then P0:=(OrderC(P0) div n)*P0; end if; 
x0:=P0[1];
S:={x0};
for i:= 1 to (n-1) do S:=S join {(i*P0)[1]}; end for;
test:=false;
i:=1;
while test eq false do if lst[i] in S then i:=i+1;  else test:=true; z0:=lst[i]; end if; end while;
z0:=K!z0;
E1:=ChangeRing(E,K);
P1:=Points(E1, z0)[1];
P0:=E1!P0;
return K,P0,P1;
end function;


//The following function tests whether an polynomial has any factors of degree dividing d
// over F_q, where q is coprime to n
// truth =  the polynomial does not have a factor of degree dividing d
// one has to also input the defining field.

function poltest(p,n,d)  
// local truth,truth2,i,fac,j,b,lc,tmp,deg,f,RF,coef,zz,pg,K;
//print "poltest entered";
truth:=false;
K:=FieldOfFractions(BaseRing(p));
i:=1;
coeff:=Coefficients(p);
b:=[Denominator(t):t in coeff];
lc:=LCM(b);
//print "computing ring of integers";
O:=RingOfIntegers(K);
//print "computed ring of integers";
deg:=Degree(AbsoluteField(K));
coef:=Coefficients(p);
while (i le 50) and (truth eq false) do
	q:=NthPrime(i);
	fac:=Factorization(q*O);
	tmp:=deg;
	tmppr:=fac[1][1];
	for j:=1 to #fac do 
		if Degree(fac[j][1]) lt tmp then 
			tmppr:=fac[j,1]; tmp:=Degree(fac[j][1]); 
		end if; 
	end for;
	if IsCoprime(q,n) and IsCoprime(q,lc) then
		RF,f:=ResidueClassField(tmppr);
		zz:=[f(bd):bd in coef];
		p1:=Polynomial(zz);
		fac:=Factorization(p1);
		truth2:=false;
		for j:=1 to #fac do
			if d mod Degree(fac[j][1]) eq 0 then truth2:=true; end if; 
		end for;
		if truth2 eq false then truth:=true; end if;
	end if;
	i:=i+1;
end while;
//print "poltest exited";
return truth;
end function;




function DivisionFieldHard(E,n)
// E an elliptic curve over a number field K
// n an odd integer
// we assume that there is no point of order n on E defined over K, so we have to add x(P_1),y(P_1) and x(P_2) and \zeta_n and 
// all the x-coordinates of the n-division polynomaial (by Bandini-Paladino, JNT, Theorem 1.1.)
// This isn't implemented properly yet
// local K,i,p,fac,lst,K1,x0,E2,pts,a1,a2,a3,a4,a6,y2,pts,P0,test,S,z0,P1;
p:=PrimitiveDivisionPolynomial(E,n);
fac:=Factorization(p);
//pos:=false;
i:=0;
K:=BaseField(E);
K1,lst:=SplittingField(p);
_<x>:=PolynomialRing(K1);
fac2:=Factorization(x^n-1);
K2:=ext<K1|fac2[#fac2][1]>;
x0:=K2!lst[1];
E2:=ChangeRing(E,K2);
pts:=Points(E2,x0);
if #pts ge 1 then P0:=pts[1]; K:=AbsoluteField(K2); E2:=ChangeRing(E,K); P0:=E2!P0; else 
	a1,a2,a3,a4,a6:=Explode(aInvariants(E2));
	y2:=Polynomial([-(x0^3+a2*x0^2+a4*x0+a6),(a1*x0+a3),1]);
	booIrre:=IsIrreducible(y2); 
	assert booIrre;
	K:=AbsoluteField(ext<K2|y2>);
	x0:=K!lst[1];
	E2:=ChangeRing(E,K);
	pts:=Points(E2,x0);
	P0:=pts[1];
end if;
test:=false;
i:=1;
S:={x0};
while test eq false do if lst[i] in S then i:=i+1;  else test:=true; z0:=lst[i]; end if; end while;
z0:=K!z0;
P1:=Points(E2, z0)[1];
return K,P0,P1;
end function;

// The following function returs Phi_n -see Washington's book.

function PhiDivisionPolynomial(E,n)
// local K,p,p1,p2,g1,g2,t1,t2,g,t;
K<x>:=PolynomialRing(BaseField(E));
p,g,t:=DivisionPolynomial(E,n);
p1,g1,t1:=DivisionPolynomial(E,n-1);
p2,g2,t2:=DivisionPolynomial(E,n+1);
if IsOdd(n) then
	return x*p^2-g1*g2*t1;
else
	return x*t*g^2-p1*p2;
end if;
end function;



//this function  takes 2 elements and determines |<P,Q>|

function nmreltsgen(P,Q)
// local i,j,ordP,ordQ,set;
set:={};
ordP:=OrderC(P);
ordQ:=OrderC(Q);
for i:=1 to ordP do
	for j:=1 to ordQ do
		set:=set join {i*P+j*Q};
	end for;
end for;
return #set;
end function;



// We suppose that as input we receive the set R of primes ell there will be ell-primary torsion growth
// The output will be a list of (K,T,gen,div), where K is a number field with primitve torsion growth, and T is the appropriate torsion group, gen is the set that generates the p-power torsion subgroup, and div is true if K is a ell^n-division field for some n
// Ge is assumed to be the order of the mod l represenation
// F will be the set of Q(E[p^n]) and A will be the set of all the other fields. 

function step2(E,ell,d,Ge)
// local A,F,T,m,P0,E1,test,K,K1,P,Q, p, fac, j, a, pts, divcomputed, d1, deg, cntr,startnmr, new, primdiv, additioncheck,tr,mapf, alp,alp2,L, Pt, ptsnew, checkiso, varbreak, temppair, ordp, ordq, B, C,sPT, PT, Knew, Knew2;
A:=<>;
F:=<>;
jI:=jInvariant(E);
if (ell eq 37) and d lt 36 then
	if (d mod 12 eq 0) and (jI eq -9317) then
		// tr,dd:=IsQuadraticTwist(E, EllipticCurve([1,1,1,-8,6])); ENRIQUE
		tr,dd:=IsQuadraticTwist(E, EllipticCurve([-1155,16450]));
		if tr then
			_<x>:=PolynomialRing(Rationals());
			//L<alp>:=NumberField(x^6-85*x^5+435*x^4-750*x^3+400*x^2+125*x-125); ENRIQUE
			L<alp>:=NumberField(x^6 - 210*x^5 - 8085*x^4 + 125300*x^3 + 4251975*x^2 - 16133250*x - 408849875);
			_<y>:=PolynomialRing(L);
			p2:=dd*y^2-(alp^3-1155*alp + 16450);
			// in this case we will not compute the point of order 37, as it won't be relevant in the program any more
			A:=<<AbsoluteField(ext<L|p2>), <37>, <>>>;
		end if;
	end if;
else
	if (ell eq 17 and (d lt 288) ) and ((d mod 8 ne 0) or ((jI ne -17*373^3/2^17) and (jI ne (-17^2*101^3)/2))) then
		A:=<>;
	else 
		if ell ge 3 then		 
			T,m:=pPowerTorsion(E,ell);
			if (#T mod ell) eq 0 then 
				Ptemp:=m(T.1);
				mult:= OrderC(Ptemp) div ell;
				P0:=mult*Ptemp;
				assert OrderC(P0) eq ell;
				A:=A cat <<RationalsAsNumberField(), <ell,1>, <P0>>>;
				if (d mod Ge eq 0) then
					K1,P,Q:=DivisionFieldEasy(E,ell);
					A:=A cat <<K1, <ell,ell>, <P,Q>>>;
				end if;
			else
				lst:=<b: b in PolDegree0(E, ell,1, d)>;
				for i:=1 to #lst do 
					if (Degree(lst[i,2]) ge 2) then 
						K:=lst[i,2];
						E1:=ChangeRing(E,K);
						//p:=PrimitiveDivisionPolynomial(E1,ell);
						//if Degree(K) eq 1 then  
							//fac:=Factorization(p:DegreeLimit:=d);
						//else 
						//	fac:=Factorization(p);
						//end if;
						//assert Degree(fac[1][1]) eq 1;
						//test:=true;
						j:=1;
						a:=lst[i,1];
						pts:=Points(E1,a);
						assert #pts ge 1;
						P0:=pts[1];
						//while test and (j le #fac) do 
							//print "while3";
						//	if Degree(fac[j][1]) eq 1 then
						//		a:=Roots(fac[j][1])[1,1];
						//		pts:=Points(E1,a);
						//		if #pts ge 1 then
						//			test:=false; 
						//			P0:=pts[1];
						//		end if;
						//	end if;
						//	j:=j+1;
						//end while;
						//assert test eq false;
						A:=A cat <<K,<ell,1>,<P0> > >; 
					end if;
				end for;
				divcomputed:=false;
				for i:= 1 to #A do
					if Degree(A[i][1]) eq Ge then 
						K1,P,Q:=DivisionFieldHard(E,ell);
						E2:=ChangeRing(E,K1);
						A:=replace(A,i, <K1,<ell,ell>,<E2!P,E2!Q>>);
						divcomputed:=true;
					end if;
				end for;
				if (divcomputed eq false) then
					if (d mod Ge) eq 0 then
						K,P0,Q0:=DivisionFieldHard(E, ell);
						A:=A cat <<K,<ell,ell>,<P0,Q0>>>;
						F:=F cat <<K,<ell,ell>,<P0,Q0>>>;
					end if;
				end if;
			end if;
		else 
			assert ell eq 2;
			T,m:=TwoTorsionSubgroup(E);
			if #T eq 2 then
				A:=A cat < <RationalsAsNumberField(), <2,1>, <m(T.1)>> >;
				if d mod 2 eq 0 then
					K1,P,Q:=DivisionFieldEasy(E,ell);
					A:=A cat <<K1, <ell,ell>, <P,Q>>>;
					F:=F cat <<K1, <ell,ell>, <P,Q>>>;
				end if;		
			else 
				if #T eq 4 then
					A:=A cat < <RationalsAsNumberField(), <2,2>, <m(T.1), m(T.2)>> >;
					F:=F cat < <RationalsAsNumberField(), <2,2>, <m(T.1), m(T.2)> > >;
				else 
					assert #T eq 1;
					if d mod 3 eq 0 then
						if Ge eq 6 then
							ptmp:=DivisionPolynomial(E,2);
							K<w>:=NumberField(ptmp);
							E1:=ChangeRing(E,K);
							pts:=Points(E1,w);
							assert #pts eq 1;
							P:=pts[1];
							A:= A cat  < <K, <2,1>, <P>> >;
							if d mod 6 eq 0 then
							K,P,Q:=DivisionFieldHard(E,2);
							A:= A cat < <K, <2,2>, <P,Q>> >;
							F:= F cat < <K, <2,2>, <P,Q>> >;
							end if;	
						else
							assert Ge eq 3;
							K,P,Q:=DivisionFieldHard(E,2);
							A:= A cat < <K, <2,2>, <P,Q>> >;
							F:= F cat < <K, <2,2>, <P,Q>> >;
						end if;
					end if;
				end if;
			end if;
		end if;
	end if;
end if;
// important note: points are always stored such that the order of the second is not larger than the order of the first

k:=2;
new:=true;
while new and (ell ne 5 or k le 2 or d ge 50) and (ell lt 11 or d ge 55) and (ell ne 7 or d ge 42) do
	if booprint then printf "k=%o\n",k; end if;
	startnmr:=#A;
	new:=false;
	primdiv:=PrimitiveDivisionPolynomial(E,ell^k);
	if poltest(primdiv,ell,d) eq false then
		if booprint then printf "test passed %o\n", k; end if;
		nmr:=#A;
		for i:=1 to nmr do
			if booprint then printf "i=%o\n",i; end if;
			K:=A[i][1];
			deg:=Degree(K);
			assert (d mod deg eq 0);
			E1:=ChangeRing(E,K);
			temppair:=A[i][3];
			d1:=Numerator(d/deg);
			if #temppair eq 1 then
				if booprint then print "entering #temppair eq 1"; end if;
				P:=temppair[1];
				if OrderC(P) eq ell^(k-1) then
					if ell ge 3 then
						poly:=PolynomialRing(K)!(PhiDivisionPolynomial(E1,ell)-P[1]*DivisionPolynomial(E1,ell)^2);
					else
						assert ell eq 2;
						poly:=PolynomialRing(K)!(PhiDivisionPolynomial(E1,ell)-P[1]*DivisionPolynomial(E1,ell));
					end if;
					testingvariable:=false;
					if booprint then printf "%o", "about to do poltest..."; t0 := Cputime(); end if;
					if ell gt 2 and Degree(K) eq 1 then testingvariable:=poltest(poly,ell,d1); end if;
					if booprint then printf "poltest passed in %o seconds\n", Cputime(t0); end if;
					if testingvariable eq false then
						fac:=Factorization(poly: DegreeLimit := d1);
						if #fac ge 1 then
							for l:=1 to #fac do
								current:=fac[l,1];
								if d1 mod Degree(current) eq 0 then
									if Degree(current) ge 2 then
										Knew2<alp2>:=ext<K|current>;
										Knew:=AbsoluteField(Knew2);
										alp:=Knew!alp2;
										//polynom:=AbsoluteMinimalPolynomial(alp2);
										//alp:=Roots(PolynomialRing(Knew)!polynom)[1,1];
										//print "about to compute isomorphism";
										//tr,mapf:=IsIsomorphic(Knew2,Knew);
										//print "done";
										//alp3:=mapf(alp2);
										//alp3,alp;
										//assert tr;
									else
										Knew:=K;
										alp:=Roots(current)[1,1];
									end if;
									checkiso:=false;
									varbreak:=true;
									if ell^k eq 4 then
										E2:=ChangeRing(E1,Knew);
										L,Pt:=NFext(E2,alp);
										if (Degree(L) eq 2*Degree(Knew)) and IsOdd(d1) or d mod Degree(L) ne 0 then
											varbreak:=false;
										else
										Knew:=L;
										alp:=Pt[1];
										end if;
									end if;
									if varbreak then
										new:=true;
										checkiso:=false;
										if Degree(Knew) gt Degree(K) then
											for st in A do
												if IsIsomorphicNF(st[1],Knew) then 
													if st[2,1] eq ell^k then 
														checkiso:=true;
													end if;
												end if; 
											end for;
										else
											checkiso:=true;
											assert IsIsomorphicNF(K,Knew);
											E2:=ChangeRing(E1,Knew);
											ptsnew:=Points(E2,alp);
											assert #ptsnew ge 1;
											Pnew:=ptsnew[1];
											A[i][2]:=<ell^k,1>; A[i][3]:=<E2!Pnew>;
											assert OrderC(E2!Pnew) eq ell^k;
										end if;
										if checkiso eq false then 
											Krep,frep:=SwitchToOptimized(Knew);
											E2:=ChangeRing(E1,Krep);
											ptsnew:=Points(E2,frep(alp));
											Knew:=Krep;
											assert #ptsnew ge 1;
											Pnew:=ptsnew[1];
											assert OrderC(E2!Pnew) eq ell^k;
											A:=A cat < <Krep, <ell^k,1>,<E2!Pnew>> >;
											// now we need to check whether we picked up any other torsion while picking up a point of order p^k 
											if #F ge 1 then
												cntr:=1;
												repeat
													additioncheck:=true;
													Fdiv:=F[cntr,1];
													if (Degree(Knew) mod Degree(Fdiv) eq 0) and (IsCoprime(Degree(Knew) div Degree(K), Degree(Fdiv)) eq false) then
														
														if IsSubfieldNF(Fdiv, Knew) then 
															if booprint then printf "%o","about to compute is subfield (3).."; t0 := Cputime(); end if;
															tr,mapf:=IsSubfieldNF(Fdiv, Knew);
																if booprint then printf "done in %o seconds\n",Cputime(t0); end if;
														else 
															tr:=false;
														end if;
														
														if tr then
															additioncheck:=false;
															P1:=E2![mapf(F[cntr,3][1][1]), mapf(F[cntr,3][1][2])];
															if nmreltsgen(E2!Pnew,E2!P1) eq ell^(k+cntr) then
																A:=replace(A,#A, <A[#A,1],<ell^k,ell^cntr>,<E2!Pnew,E2!P1>>);
																assert OrderC(E2!Pnew) eq ell^k;
																assert OrderC(E2!P1) eq ell^cntr;
																//A[#A,2]:=<ell^k,ell^cntr>; A[#A,3]:=<E2!Pnew,E2!P1>;
															else 
																P2:=E2![mapf(F[cntr,3][2][1]), mapf(F[cntr,3][2][2])];
																assert nmreltsgen(E2!Pnew,E2!P2) eq ell^(k+cntr);
																assert IsIsomorphicNF(Knew,A[#A,1]);
																A:=replace(A,#A, <A[#A,1],<ell^k,ell^cntr>,<E2!Pnew,E2!P2>>);
																assert OrderC(E2!Pnew) eq ell^k;
																assert OrderC(E2!P2) eq ell^cntr;
																//A[#A,2]:=<ell^k,ell^cntr>; A[#A,3]:=<E2!Pnew,E2!P2>;
															end if;
														end if;
													end if;
													cntr:=cntr+1;
												until additioncheck or (cntr ge k) or (cntr gt #F);
											end if;
										end if;
									end if;
								end if;
							end for;
						end if;
					end if;
				end if;
			else
				if booprint then print "entering #temppair eq 2"; end if;
				assert #temppair eq 2;
				P:=temppair[1];
				Q:=temppair[2];
				ordp:=OrderC(P);
				ordq:=OrderC(Q);
				assert ordp-ordq ge 0;
				//nmreltsgen(P,Q), A[i,2,1], A[i,2,2], ordp, ordq;
				assert ordp eq A[i,2,1];
				assert nmreltsgen(P,Q) eq A[i,2,1]*A[i,2,2];				
				//print "asserts passed";
				if (ordp eq ell^(k-1)) then 
					setofpts:={};
					for num:=0 to (ordq-1) do setofpts:=setofpts join {<P+num*Q,num>}; 
						if OrderC(P+num*Q) ne ell^(k-1) then print "greska:",i,nmreltsgen(P,Q),ordp,A[i,2,1], ell^k, num, OrderC(P+num*Q), ell^(k-1), nmreltsgen(P,Q); end if;
					end for; 
					if ordp eq ordq then
						setofpts:=setofpts join {<Q,ordq+1>};
					end if;
					for sPT in setofpts do
						if A[i,2,1] lt ell^k then
							PT:=sPT[1];
							if ell ge 3 then
								poly:=PolynomialRing(K)!(PhiDivisionPolynomial(E1,ell)-PT[1]*DivisionPolynomial(E1,ell)^2);
							else
								assert ell eq 2;
								poly:=PolynomialRing(K)!(PhiDivisionPolynomial(E1,ell)-PT[1]*DivisionPolynomial(E1,ell));
							end if;
							testingvariable:=false;
							if ell gt 2 and Degree(K) eq 1 then testingvariable:=poltest(poly,ell,d1); end if;
							if testingvariable eq false then
								fac:=Factorization(poly: DegreeLimit := d1);
								if #fac ge 1 then
									for l:=1 to #fac do
									current:=fac[l,1];
									if d1 mod Degree(current) eq 0 then
										if Degree(current) ge 2 then
											Knew2<alp2>:=ext<K|current>;
											Knew:=AbsoluteField(Knew2);
											alp:=Knew!alp2;
											//polynom:=AbsoluteMinimalPolynomial(alp2);
											//alp:=Roots(PolynomialRing(Knew)!polynom)[1,1];
											//print "about to compute isomorphism";
											//tr,mapf:=IsIsomorphic(Knew2,Knew);
											//alp3:=mapf(alp2);
											//alp3,alp;
											//print "done";
											//alp:=mapf(alp2);
											//assert tr;
										else
											Knew:=K;
											alp:=Roots(current)[1,1];
										end if;
										checkiso:=false;
										varbreak:=true;
										if ell^k eq 4 then
											E2:=ChangeRing(E1,Knew);
											L,Pt:=NFext(E2,alp);
											//Order(Pt),ell^k;
											if (Degree(L) eq 2*Degree(Knew)) and IsOdd(d1) or d mod Degree(L) ne 0 then
												varbreak:=false;
											else
												if Degree(L) eq 2*Degree(Knew) then
													Kold:=Knew;
													Knew:=L;
													//print "computing subfield";
													//tr2,mapf2:=IsSubfieldNF(Kold, Knew);
													//print "computed subfield";
													//assert tr2;
													alpold:=alp;
													//alp:=mapf2(alpold);
													alp:=Pt[1];
												else
													assert Degree(L) eq Degree(Knew);
												end if;
											end if;
											
										end if;
										if varbreak then
											new:=true;
											if Degree(Knew) gt Degree(K) then
												for st in A do
													//print "checking is isomorphic";
													if IsIsomorphicNF(st[1],Knew) then 
														if st[2,1] ge ell^k and (st[2,2] ge A[i,2,2]) then 
															checkiso:=true;
														end if;
													end if; 
													//print "checked isomorphic";
												end for;
											else
												checkiso:=true;
												//print "checking is isomorphic";
												assert IsIsomorphicNF(K,Knew);
												//print "checked isomorphic";
												E2:=ChangeRing(E1,Knew);
												ptsnew:=Points(E2,Knew!alp);
												assert #ptsnew ge 1;
												Pnew:=ptsnew[1];
												//Order(Pnew),OrderC(Pnew);
												assert OrderC(Pnew) eq ell^k;
												A[i][2]:=<ell^k,A[i][2][2]>; 
												Qtmp:=A[i,3,2];
												Ptmp:=A[i,3,1];
												if nmreltsgen(E2!Pnew,E2!Qtmp) eq A[i,2,1]*A[i,2,2] then
													A[i][3]:=<E2!Pnew, E2!Qtmp>;
													//k,A[i,2,2],OrderC(E2!Qtmp), A[i,2,1], OrderC(Pnew);
													assert A[i,2,2] eq OrderC(E2!Qtmp) and A[i,2,1] eq OrderC(E2!Pnew);
													assert A[i,2,1]*A[i,2,2] eq nmreltsgen(E2!Pnew,E2!Qtmp);

												else 
													if (nmreltsgen(E2!Pnew,E2!Qtmp) gt A[i,2,1]*A[i,2,2]) and (nmreltsgen(E2!Pnew,E2!Ptmp) gt A[i,2,1]*A[i,2,2]) then
														A[i,2,2]:=ell^k;
														A[i][3]:=<E2!Pnew, E2!Ptmp>;
														assert A[i,2,2] eq OrderC(E2!Qtmp) and A[i,2,1] eq OrderC(Pnew);
														assert A[i,2,1]*A[i,2,2] eq nmreltsgen(E2!Pnew,E2!Ptmp);
													else
														assert A[i][2][2] eq ell^(k-1);
														A[i][3]:=<E2!Pnew, E2!Ptmp>;
													//print k,A[i,2,1],A[i,2,2],OrderC(E2!Pnew), OrderC(E2!Ptmp);
														assert A[i,2,1]*ell^(k-1) eq nmreltsgen(E2!Pnew,E2!Ptmp);
														assert A[i,2,1]*A[i,2,2] eq nmreltsgen(E2!Pnew,E2!Ptmp);
													end if;
												end if;
											end if;
											if checkiso eq false then 
												Krep,frep:=SwitchToOptimized(Knew);
												Knew:=Krep;
												E2:=ChangeRing(E1,Krep);
												ptsnew:=Points(E2,frep(alp));
												//E2:=ChangeRing(E1,Knew);
												//ptsnew:=Points(E2,alp);
												assert #ptsnew ge 1;
												Pnew:=ptsnew[1];
												//Order(Pnew), ell^k;
												assert OrderC(E2!Pnew) eq ell^k;
												num:=sPT[2];
												Prep:=E2![frep(P[1]),frep(P[2])];
												Qrep:=E2![frep(Q[1]),frep(Q[2])];
												if (num eq ordq+1) then 
													A:=A cat < <Krep, <ell^k,ell^(k-1)>,<E2!Pnew,E2!Prep>> >;
													assert A[#A,2,2] eq OrderC(E2!Prep) and A[#A,2,1] eq OrderC(E2!Pnew);
													assert A[#A,2,1]*A[#A,2,2] eq nmreltsgen(E2!Pnew,E2!Prep); 
												else
													A:=A cat < <Krep, <ell^k,ordq>,<E2!Pnew,E2!Qrep>> >;
													assert A[#A,2,2] eq OrderC(E2!Qrep) and A[#A,2,1] eq OrderC(E2!Pnew);
													//A[#A,2,1],A[#A,2,2], OrderC(E2!Pnew),OrderC(E2!Qrep),nmreltsgen(E2!Pnew,E2!Qrep);  
													assert A[#A,2,1]*A[#A,2,2] eq nmreltsgen(E2!Pnew,E2!Qrep);
												end if;
												//now we check whether we've picked up any new torsion
												cntr:=Round(Log(ell,A[#A,2][2]))+1;
												additioncheck:=false;
												while (additioncheck eq false) and (cntr le #F) do
													//print "while1";
													additioncheck:=true;
													Fdiv:=F[cntr,1];
													if (Degree(Knew) mod Degree(Fdiv) eq 0) and (IsCoprime(Degree(Knew) div Degree(K), Degree(Fdiv)) eq false) then
														
														if IsSubfieldNF(Fdiv, Knew) then 
															if booprint then printf "%o","compute subfield (2)..."; t0 := Cputime(); end if; 
															tr,mapf:=IsSubfieldNF(Fdiv, Knew);
															if booprint then printf "done in %o seconds\n",Cputime(t0); end if;
														else 
															tr:=false;
														end if;
														//tr,mapf:=IsSubfieldNF(Fdiv, Knew);
														
														if tr and A[#A,2,2] lt ell^cntr then
															additioncheck:=false;
															P1:=E2![mapf(F[cntr,3][1][1]), mapf(F[cntr,3][1][2])];
															if nmreltsgen(E2!Pnew,E2!P1) eq ell^(k+cntr) then
																A:=replace(A, #A, <A[#A,1], <ell^k,ell^cntr>, <E2!Pnew,E2!P1>>);
																//A[#A,2]:=<ell^k,ell^cntr>; A[#A,3]:=<E2!Pnew,E2!P1>;
															else 
																P2:=E2![mapf(F[cntr,3][2][1]), mapf(F[cntr,3][2][2])];
																if nmreltsgen(E2!Pnew,E2!P2) eq ell^(k+cntr) then
																	A:=replace(A, #A, <A[#A,1], <ell^k,ell^cntr>, <E2!Pnew,E2!P2>>);
																else 
																	//nmreltsgen(E2!Pnew,E2!P2),nmreltsgen(E2!Pnew,E2!P2),ell^(k+cntr);
																	F;
																end if;
																//A:=replace(A, #A, <A[#A,1], <ell^k,ell^cntr>, <E2!Pnew,E2!P2>>);
																//A[#A,2]:=<ell^k,ell^cntr>; A[#A,3]:=<E2!Pnew,E2!P2>;
															end if;
														end if;
													end if;
													cntr:=cntr+1;
												end while; 
											end if;
										end if;
									end if;
									end for;
								end if;
							end if;
						end if;
					end for;
				end if;
			end if;
		end for;
	end if;
	// we've now computed all the fields with points of order ell^k; it reamins to check whether the ell^k-division field is of degree dividing d 
	// we do this by only checking fields that have torsion (ell^k, ell^(k-1))
	
	
	if k gt #F then
		//print "entering division field part";
		for c:=1 to #A do
			//c, A[c];
			if A[c,2,1] eq ell^k and A[c,2,2] eq ell^(k-1) and k gt #F then
				PT:=A[c,3,2];
				E1:=ChangeRing(E,A[c,1]);
				if ell ge 3 then
					poly:=PolynomialRing(A[c,1])!(PhiDivisionPolynomial(E1,ell)-PT[1]*DivisionPolynomial(E1,ell)^2);
				else
					assert ell eq 2;
					poly:=PolynomialRing(A[c,1])!(PhiDivisionPolynomial(E1,ell)-PT[1]*DivisionPolynomial(E1,ell));
				end if;
				//print "entered 1 if and passed 1 if";
				d1:= d div Degree(A[c,1]);
				testingvariable:=false;
				if ell gt 2 and Degree(K) eq 1 then testingvariable:=poltest(poly,ell,d1); end if;
				if testingvariable eq false then
					//print "factoring division polynomial";
					fac:=Factorization(poly: DegreeLimit := d1);
					//print "factored division polynomial";
					if #fac ge 1 then
						K:=A[c,1];
						for l:=1 to #fac do
							if k gt #F then 
								current:=fac[l,1];
								if d1 mod Degree(current) eq 0 then
									new:=true;
									varbreak:=true;
									if ell^k eq 4 then
										if Degree(current) gt 1 then
											Knew2<alp2>:=ext<K|current>;
											Knew:=AbsoluteField(Knew2);
											alp:=Knew!alp2;
											//polynom:=AbsoluteMinimalPolynomial(alp2);
											//alp:=Roots(PolynomialRing(Knew)!polynom)[1,1];
											//print "testing isomorphism";
											//tr,mapf:=IsIsomorphic(Knew2,Knew);
											//alp3:=mapf(alp2);
											//alp3,alp;
											//print "done";
											//assert tr;
											//alp:=mapf(alp2);
										else
											alp:=Roots(current)[1,1];
											Knew:=K;											
										end if;
										E2:=ChangeRing(E,Knew);
										L,Pt:=NFext(E2,alp);
										Knew:=L;
										E2:=ChangeRing(E,Knew);
										alp:=(E2!Pt)[1];
										if (Degree(L) eq 2*Degree(Knew)) and IsOdd(d1) or d mod Degree(L) ne 0 then
											varbreak:=false;
										else
											Knew:=L;
											alp:=Pt[1];
										end if;
									else 
										if Degree(current) gt 1 then 
											Knew2<alp2>:=ext<K|current>;
											Knew:=AbsoluteField(Knew2);
											alp:=Knew!alp2;
											//polynom:=AbsoluteMinimalPolynomial(alp2);
											//alp:=Roots(PolynomialRing(Knew)!polynom)[1,1];
											//print "about to compute isomorphism";
											//tr,mapf:=IsIsomorphic(Knew2,Knew);
											//alp3:=mapf(alp2);
											//alp3,alp;
											//print "done";
											//assert tr;
											//alp:=mapf(alp2);
											//assert Evaluate(poly,alp2) eq 0;
										else
											Knew:=K;
											alp:=Roots(current)[1,1];
										end if;
											E2:=ChangeRing(E,Knew);
										
									end if;	
									if varbreak then
										Krep,frep:=SwitchToOptimized(Knew);
										E2:=ChangeRing(E1,Krep);
										ptsnew:=Points(E2,frep(alp));
										
										//ptsnew:=Points(E2,alp);
										assert #ptsnew ge 1;
										Pnew:=ptsnew[1];
										assert OrderC(Pnew) eq ell^k;
										assert OrderC(A[c,3,1]) eq ell^k;
										assert nmreltsgen(E2!Pnew,E2!A[c,3,1]) eq ell^(2*k);
										assert(Degree(Knew) le d);
										assert k gt #F;
										A:=A cat < <Krep, <ell^k,ell^k>,<E2![frep(A[c,3,1,1]),frep(A[c,3,1,2])],E2!Pnew>> >;
										F:=F cat < <Krep, <ell^k,ell^k>,<E2![frep(A[c,3,1,1]),frep(A[c,3,1,2])],E2!Pnew>> >;
									end if;
								end if;
							end if;
						end for;					
					end if;
				end if;
			end if;
		end for;
		//print "exiting division field";
	end if;
	k:=k+1;
end while;
B:=<>;
for i:=1 to #A do
	check:=true;
	for j:= 1 to #B do
	if IsIsomorphicNF(A[i,1], B[j,1]) then
		if A[i,2,1] ge B[j,2,1] and A[i,2,2] ge B[j,2,2] 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;
			check:=false;
		else
			check:=false;		
		end if;
	end if;
	end for;
	if check then 
		B:= B cat <A[i]>;
	end if;
end for;
return B;

end function;

// todo Optimized representation? load "program.m"; E:=EllipticCurve("101388b1"); TorsionGrowth(E,8);
