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


////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Andrew Sutherland "Computing images of Galois representations attached to elliptic curves", Forum of Mathematics, Sigma 4 (2016) e4 (79 pages), http://dx.doi.org/10.1017/fms.2015.33
// https://math.mit.edu/~drew/galrep/
//
load "galrep_Sutherland/galrep.m";
load "galrep_Sutherland/subgroups.m";
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// booprintGal:=false;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

function IsSubfieldNF(K1,K2)
// K1 and K2 number fields
// the function return if K1 is a subfield  of K2 or not
  if not IsDivisibleBy(Degree(K2),Degree(K1)) then return false; end if;
  p1:=DefiningPolynomial(K1);
  p2:=DefiningPolynomial(K2);
  if Degree(K1) eq 1 and Degree(K2) eq 1 then return true, IdentityHomomorphism(RationalsAsNumberField()); end if;
  raiz:=Roots(ChangeRing(p1,K2));
  if #raiz eq 0 then return false,0; end if;
  s:=raiz[1,1];
  assert Evaluate(ChangeRing(p1,K2),s) eq 0;
  return true, hom<K1->K2 | s>; 
end function;

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

function IsIsomorphicNF(K1,K2)
// K1 and K2 number fields
// the function return if K1 is isomorphic to K2 or not
  if Degree(K1) ne Degree(K2) then return false,0; end if;
  boo,map:=IsSubfieldNF(K1,K2);
  return boo,map;
end function;


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

function CompositumNF(K1,K2)
// K1 and K2 absolute number fields
// returns K1K2 (as an absolute number field)
p2:=(PolynomialRing(K1))!(DefiningPolynomial(K2));
fac:=Factorization(p2);
return AbsoluteField(ext<K1|fac[1][1]>);
end function;

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

function aux_y(E,a,d)
// E an elliptic curve over Q
// a an algebraic number
// d a positive integer
//
// OUTPUT:
// Defining polynomials of the number field Q(a,b) such that (a,b) in E *ONLY* in the case that [Q(a,b):Q] divides d.
// If E:y^2=x^3+A*x+B then b^2=a^3+A*a+b

  K:=FieldOfFractions(Parent(a));
  if not IsDivisibleBy(d,Degree(K)) then return false,0,0; end if;
  if Degree(K) eq 1 then K:=RationalsAsNumberField(); end if;
  a1,a2,a3,a4,a6:=Explode(aInvariants(BaseChange(E,K)));
  y2:=Polynomial([-(a^3+a2*a^2+a4*a+a6),(a1*a+a3),1]);
  booIrre:=IsIrreducible(y2);
  //if (IsOdd(d) and booIrre) or (booIrre and Degree(K) eq d) or (not IsDivisibleBy(d,Degree(K))) then return false,0,0; end if;
  if not booIrre then
    pts:=Points(ChangeRing(E,K),K!a);
	assert #pts ge 1;
    return true,K,a;
   else
    L:=AbsoluteField(NumberField(y2));
    if not IsDivisibleBy(d,Degree(L)) then return false,0,0; end if; 
    pts:=Points(ChangeRing(E,L),L!a);
	assert #pts ge 1;
    return true,L,L!a; 
  end if;
end function;


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


function PolDegree0(E,p,n,d)
// Input: E an elliptic curve over Q
// 		  p a prime
//        n and d positive integers
// OUTPUT:  Tuple of pairs <alpha,K> such that 
//          P=(alpha,beta) is a point of order p^n
//          and K=Q(P)
//  local a,fpn,Lpnd,fa;
  fpn:=PrimitiveDivisionPolynomial(E,p^n);   
  Lpnd:=[f[1] : f in Factorization(fpn : DegreeLimit := d) | IsDivisibleBy(d,Degree(f[1]))];

  L:=<>;
  for f in Lpnd do
    if Degree(f) eq 1 then
      a:=Roots(f)[1,1]; 
     else
      K<a>:=NumberField(f);
    end if;
    boo,La,aa:=aux_y(E,a,d);
    if boo then
        L:=Append(L,<aa,La>);
    end if;
  end for;
  if #L eq 0 then return <>; end if;
  
  LL:=<L[1]>;
  for k in [2..#L] do
    if {IsIsomorphicNF(L[k][2],K[2]) : K in LL} eq {false} then
      LL:=Append(LL,L[k]); 
    end if; 
  end for;
  return LL;
end function;

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

function aux_IsSurjective(E,p)
// Input: E an elliptic curve without CM over Q and p>11 a prime
// Output: If return true then G_E(p)=GL_2(Fp).
//		   If return false then G_E(p)\ne GL_2(Fp)
	assert p ge 13;
	assert not HasComplexMultiplication(E);
	jE:=jInvariant(E);
	assert not [p,jE] in {[17,-17*373^3/2^17],[17,-17^2*101^3/2],[37,-7*11^3],[37,-7*137^3*2083^3]};
// Theorem 1.5
// D. Zywina, "On the surjectivity of mod ell representations associated to elliptic curves". arXiv: 1508.07661
	if p gt 13 then
		den:=Denominator(jE);
		L:=Factorization(den);
		if not &and[ IsDivisibleBy(f[2],p) : f in L] or &and[IsDivisibleBy(q[1]^2-1,p) : q in L] then return true; end if;
	end if;

// Proposition 6.4
// D. Zywina, "On the surjectivity of mod ell representations associated to elliptic curves". arXiv: 1508.07661
	bound:=10;
	c:=0; 
	q:=2;
	NE:=Conductor(E);
	while c lt bound do
		q:=NextPrime(q);
		if not IsDivisibleBy(p*NE,q) then
			aq:=TraceOfFrobenius(E,q);
			bq:=aq^2-4*q;
			if (not IsDivisibleBy(aq,p))and (not IsDivisibleBy(bq,p)) and LegendreSymbol(bq,p) eq 1 then
				return true;
			end if;
			c:=c+1;
		end if;
	end while;
// Proposition 6.3 (ii) 2 last items. The first one is Prop 6.4 for p=13
// D. Zywina, "On the surjectivity of mod ell representations associated to elliptic curves". arXiv: 1508.07661
	if p eq 13 then
		bound:=10;
		c:=0; 
		q:=2;
		NE:=Conductor(E);
		while c lt bound do
			q:=NextPrime(q);
			if not IsDivisibleBy(p*NE,q) then
				aq:=TraceOfFrobenius(E,q);
				if not IsDivisibleBy(aq,p) then
					return true;
				 else 
				 	d:=Integers(13)!(aq^2/p);
				 	if not d in {0,1,2,4} and d^2 -3*q+1 ne 0 then return true; end if; 
				end if;
				c:=c+1;
			end if;
		end while;	
	end if;
// Sutherland "galrep.m";
	if p eq 13 and #GaloisGroup(DivisionPolynomial(E,13)) eq 13*(13-1)*(13^2-1)/2 then return true; end if;
	GEp:=GaloisImage(E,p);
	if Order(GEp) eq p*(p-1)*(p^2-1) then return true; end if;
	print "Counterexample to Uniformity Conjecture";
	print "Uniformity Conjecture is FALSE";
// Sutherland  "subgroups.m"
	print E,p,GL2SubgroupLabel(GEp); 
	return false;
end function;


////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// D. Zywina, "On the surjectivity of mod ell representations associated to elliptic curves". arXiv: 1508.07661
//
// Appendix A. Some code:
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

ExceptionalSet:=function(E)
// Input: E an elliptic curve without CM over Q and p>11 a prime
// Output: If return a set that contains the set of exceptional primes of E (i.e. primes p such that the mod p Gal. rep is non-surjective)
	assert not HasComplexMultiplication(E);
	j:=jInvariant(E); den:=Denominator(j);
	S:={2,3,5,7,13};
	if j in {-11^2,-11*131^3} then S:=S join {11}; end if;
	if j in {-297756989/2, -882216989/131072} then S:=S join {17}; end if;
	if j in {-9317, -162677523113838677} then S:=S join {37}; end if;
	if den ne 1 then
		ispow,b,e:=IsPower(den);
		if ispow then
			P:={p: p in PrimeDivisors(e) | p ge 11};
			if P ne {} then
				g:=GCD({&*P} join {p^2-1 : p in PrimeDivisors(b)});
				S:= S join {ell : ell in PrimeDivisors(g) | ell ge 11};
			end if;
		end if;
	  else
		D:=Discriminant(E);
		Q:=PrimeDivisors( GCD(Numerator(j-1728),Numerator(D)*Denominator(D)));
		Q:=[q: q in Q | q ne 2 and IsOdd(Valuation(j-1728,q))];
		if Valuation(j,2) in {3,6,9} then Q:=[2] cat Q; end if;
		p:=2;
		alpha:=[]; beta:=[];
		repeat
		a:=0;
		while a eq 0 do
			p:=NextPrime(p); K:=KodairaSymbol(E,p);
			if K eq KodairaSymbol("I0") then
				a:=TraceOfFrobenius(E,p);
			  elif K eq KodairaSymbol("I0*") then
				a:=TraceOfFrobenius(QuadraticTwist(E,p),p);
			end if;
		end while;
		S:=S join {ell : ell in PrimeDivisors(a) | ell gt 13};
		alpha:= alpha cat [[(1-KroneckerSymbol(q,p)) div 2 : q in Q]];
		beta:= beta cat [ [(1-KroneckerSymbol(-1,p)) div 2] ];
		A:=Matrix(GF(2),alpha); b:=Matrix(GF(2),beta);
		until IsConsistent(Transpose(A),Transpose(b)) eq false;
	end if;
	return S;
end function;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  D. Zywina, "On the possible images of the mod l representations associated to elliptic curves over Q". arXiv:1508.07660
//
//  j-maps from Zywina from http://pi.math.cornell.edu/~zywina/papers/PossibleImages/
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////


F<t> := FunctionField(Rationals());
grp := recformat< level:IntegerRing(), zywina:IntegerRing(), group:GrpMat, jmap:FldFunRatUElt>;
Groups := AssociativeArray();
Groups["2Cs"] := rec<grp| level:=2, zywina:=1, group:=sub<GL(2,Integers(2)) | {[1,0,0,1]}>, jmap := 256*(t^2+t+1)^3/(t^2*(t+1)^2)>; 
Groups["2B"] := rec<grp| level:=2, zywina:=2, group:=sub<GL(2,Integers(2)) | {[1,1,0,1]}>, jmap := 256*(t+1)^3/t>;
Groups["2Cn"] := rec<grp| level:=2, zywina:=3, group:=sub<GL(2,Integers(2)) | {[1,1,1,0]}>, jmap := t^2+1728>;

Groups["3Cs"] := rec<grp| level:=3, zywina:=1, group:=sub<GL(2,Integers(3)) | {[-1,0,0,-1],[1,0,0,2]}>, jmap := 27*(t+1)^3*(t+3)^3*(t^2+3)^3/(t^3*(t^2+3*t+3)^3)>;
Groups["3Ns"] := rec<grp| level:=3, zywina:=2, group:=sub<GL(2,Integers(3)) | {[1,0,0,2],[2,0,0,1],[0,-1,1,0]}>, jmap := 27*(t+1)^3*(t-3)^3/t^3>;
Groups["3B"] := rec<grp| level:=3, zywina:=3, group:=sub<GL(2,Integers(3)) | {[1,1,0,1],[2,0,0,1],[1,0,0,2]}>, jmap := 27*(t+1)*(t+9)^3/t^3>;
Groups["3Nn"] := rec<grp| level:=3, zywina:=4, group:=sub<GL(2,Integers(3)) | {[1,0,0,-1],[-1,0,0,-1],[0,-1,1,0],[1,1,-1,1]}>, jmap := t^3>;

W3:=[ [-3*(t+1)*(t+3)*(t^2+3),-2*(t^2-3)*(t^4+6*t^3+18*t^2+18*t+9)], [], [-3*(t+1)^3*(t+9),-2*(t+1)^4*(t^2-18*t-27)],  [] ];


Groups["5Cs.4.1"] := rec<grp| level:=5, zywina:=1, group:=sub<GL(2,Integers(5)) | {[-1,0,0,-1], [1,0,0,2]}>, jmap := (t^20+228*t^15+494*t^10-228*t^5+1)^3/(t^5*(t^10-11*t^5-1)^5)>;
Groups["5Cs"] := rec<grp| level:=5, zywina:=2, group:=sub<GL(2,Integers(5)) | {[2,0,0,1], [1,0,0,2]}>, jmap := (t^2+5*t+5)^3*(t^4+5*t^2+25)^3*(t^4+5*t^3+20*t^2+25*t+25)^3/(t^5*(t^4+5*t^3+15*t^2+25*t+25)^5)>;
Groups["5Ns.2.1"] := rec<grp| level:=5, zywina:=3, group:=sub<GL(2,Integers(5)) | {[0,6,3,0],[2,0,0,2], [1,0,0,-1]}>, jmap := 5^4*t^3*(t^2+5*t+10)^3*(2*t^2+5*t+5)^3*(4*t^4+30*t^3+95*t^2+150*t+100)^3/((t^2+5*t+5)^5*(t^4+5*t^3+15*t^2+25*t+25)^5)>;
Groups["5Ns"] := rec<grp| level:=5, zywina:=4, group:=sub<GL(2,Integers(5)) | {[2,0,0,1], [1,0,0,2],[0,-1,1,0]}>, jmap := (t+5)^3*(t^2-5)^3*(t^2+5*t+10)^3/(t^2+5*t+5)^5>;
Groups["5B.4.2"] := rec<grp| level:=5, zywina:=5, group:=sub<GL(2,Integers(5)) | {[-1,0,0,-1], [1,1,0,1], [2,0,0,1]}>, jmap := ((t^4+228*t^3+494*t^2-228*t+1)^3/(t*(t^2-11*t-1)^5))>;
Groups["5B.4.1"] := rec<grp| level:=5, zywina:=6, group:=sub<GL(2,Integers(5)) | {[-1,0,0,-1], [1,1,0,1], [1,0,0,2]}>, jmap := (t^4-12*t^3+14*t^2+12*t+1)^3/(t^5*(t^2-11*t-1))>;
Groups["5Nn"] := rec<grp| level:=5, zywina:=7, group:=sub<GL(2,Integers(5)) | {[1,2,1,1],[1,6,3,1], [1,0,0,-1]}>, jmap := 5^3*(t+1)*(2*t+1)^3*(2*t^2-3*t+3)^3/(t^2+t-1)^5>;
Groups["5B"] := rec<grp| level:=5, zywina:=8, group:=sub<GL(2,Integers(5)) | {[2,0,0,1],[1,0,0,2],[1,1,0,1]}>, jmap := 5^2*(t^2+10*t+5)^3/t^5>;
Groups["5S4"] := rec<grp| level:=5, zywina:=9, group:=sub<GL(2,Integers(5)) | {[2,0,0,1], [1,0,0,2],[0,-1,1,0],[1,1,1,-1]}>, jmap :=  t^3*(t^2+5*t+40) >;

W5:=[
	[-27*(t^20+228*t^15+494*t^10-228*t^5+1), 54*(t^30-522*t^25-10005*t^20-10005*t^10+522*t^5+1)],
	[],[],[],
	[-27*(t^4+228*t^3+494*t^2-228*t+1),54*(t^6-522*t^5-10005*t^4-10005*t^2+522*t+1)],
	[-27*(t^4-12*t^3+14*t^2+12*t+1),54*(t^6-18*t^5+75*t^4+75*t^2+18*t+1)],
	[],[],[]
    	];

Groups["7Ns.3.1"] := rec<grp| level:=7, zywina:=1, group:=sub<GL(2,Integers(7)) | {[2,0,0,4], [0,2,1,0],[-1,0,0,-1]}>, jmap := F!3^3*5*7^5/2^7>;
Groups["7Ns"] := rec<grp| level:=7, zywina:=2, group:=sub<GL(2,Integers(7)) | {[3,0,0,1],[1,0,0,3],[0,-1,1,0]}>, jmap := t*(t+1)^3*(t^2-5*t+1)^3*(t^2-5*t+8)^3*(t^4-5*t^3+8*t^2-7*t+7)^3/(t^3-4*t^2+3*t+1)^7>;
Groups["7B.6.1"] := rec<grp| level:=7, zywina:=3, group:=sub<GL(2,Integers(7)) | {[-1,0,0,-1],[1,1,0,1],[1,0,0,3]}>, jmap := (t^2-t+1)^3*(t^6-11*t^5+30*t^4-15*t^3-10*t^2+5*t+1)^3/((t-1)^7*t^7*(t^3-8*t^2+5*t+1))>;
Groups["7B.6.3"] := rec<grp| level:=7, zywina:=4, group:=sub<GL(2,Integers(7)) | {[-1,0,0,-1],[1,1,0,1],[3,0,0,1]}>, jmap := (t^2-t+1)^3*(t^6+229*t^5+270*t^4-1695*t^3+1430*t^2-235*t+1)^3/((t-1)*t*(t^3-8*t^2+5*t+1)^7)>;
Groups["7B.6.2"] := rec<grp| level:=7, zywina:=5, group:=sub<GL(2,Integers(7)) | {[1,1,0,1],[3,0,0,3],[1,0,0,-1]}>, jmap := -(t^2-3*t-3)^3*(t^2-t+1)^3*(3*t^2-9*t+5)^3*(5*t^2-t-1)^3/((t^3-2*t^2-t+1)*(t^3-t^2-2*t+1)^7)>;
Groups["7Nn"] := rec<grp| level:=7, zywina:=6, group:=sub<GL(2,Integers(7)) | {[1,-4,4,1],[1,0,0,-1]}>, jmap := 64*t^3*(t^2+7)^3*(t^2-7*t+14)^3*(5*t^2-14*t-7)^3/(t^3-7*t^2+7*t+7)^7>;
Groups["7B"] := rec<grp| level:=7, zywina:=7, group:=sub<GL(2,Integers(7)) | {[3,0,0,1],[1,0,0,3],[1,1,0,1]}>, jmap := (t^2+245*t+2401)^3*(t^2+13*t+49)/t^7>;

W7:=[
    	[F!-5^3*7^3,-5^4*7^2*106],
    	[],
	[-27*(t^2-t+1)*(t^6-11*t^5+30*t^4-15*t^3-10*t^2+5*t+1),54*(t^12-18*t^11+117*t^10-354*t^9+570*t^8-486*t^7+273*t^6-222*t^5+174*t^4-46*t^3-15*t^2+6*t+1)],	
	[-27*(t^2-t+1)*(t^6+229*t^5+270*t^4-1695*t^3+1430*t^2-235*t+1),54*(t^12-522*t^11-8955*t^10+37950*t^9-70998*t^8+131562*t^7-253239*t^6+316290*t^5-218058*t^4+80090*t^3-14631*t^2+510*t+1)],
	[-27*7*(t^2-3*t-3)*(t^2-t+1)*(3*t^2-9*t+5)*(5*t^2-t-1),-54*7^2*(t^4-6*t^3+17*t^2-24*t+9)*(3*t^4-4*t^3-5*t^2-2*t-1)*(9*t^4-12*t^3-t^2+8*t-3)],
	[],
	[-27*(t^2+13*t+49)^3*(t^2+245*t+2401), 54*(t^2+13*t+49)^4*(t^4-490*t^3-21609*t^2-235298*t-823543)]
];


Groups["11B.10.4"] := rec<grp| level:=11, zywina:=1, group:=sub<GL(2,Integers(11)) | {[-1,0,0,-1],[1,1,0,1],[4,0,0,6]}>, jmap := F!-11^2>;
Groups["11B.10.5"] := rec<grp| level:=11, zywina:=2, group:=sub<GL(2,Integers(11)) | {[-1,0,0,-1],[1,1,0,1],[5,0,0,7]}>, jmap := F!-11*131^3>;

// There modular curve associated to "11Nn" is an elliptic curve over Q with positive rank so the jmap is not just a rational function in one variable. 
// Instead we use section 4.5.5 of Zywina paper
//Groups["11Nn"] :=rec<grp|level:=11, group:=sub<GL(2,Integers(11)) | {[1,-4,4,1], [1,0,0,-1]}>,jmap:=F!0 >; 
//This is a genus 1 curve and I need to think a little harder about how to work with it.
A11:=(t^5-9*t^4+17*t^3+20*t^2-73*t+43)^11;
B11:=-(t^2+3*t-6)^3*(108000*t^49 + 23793840*t^48 - 413223722*t^47 - 5377010368*t^46 + 230799738529*t^45 - 3137869050351*t^44 + 23205911712335*t^43 - 90936268647246*t^42 + 33563647471596*t^41 + 1631415220074871*t^40 - 7744726079195413*t^39 - 3218815397602111*t^38 + 236712051437217644*t^37 - 1686428698022253344*t^36 + 7984804002023063554*t^35 - 30444784135263860996*t^34 + 96849826504401032248*t^33 - 232064394883539673213*t^32 + 210175535413395353857*t^31 + 1609695806324946484826*t^30 - 11768533689837648360109*t^29 + 48291196122826259771817*t^28 - 143943931899306373170309*t^27 + 315827025781563232420857*t^26 - 421596979720485992629121*t^25 - 234929885880162547645306*t^24 + 3668241437553022801950917*t^23 - 14221091463553801024770599*t^22 + 39148264563215734730610917*t^21 - 87534472061810348609315974*t^20 + 166474240219619575379485393*t^19 - 275040771573054834247036345*t^18 + 399144725377223909937142938*t^17 - 511840960382358144595839458*t^16 + 581656165535334214665717816*t^15 - 586206578096981243980668654*t^14 + 523465655841901079370457175*t^13 - 413200824632802503354807972*t^12 + 287270832775316643952335709*t^11 - 175049577131269087795781453*t^10 + 92916572268973769104815620*t^9 - 42636417323385892254033027*t^8 + 16754292456737738144357709*t^7 - 5570911068111617263502302*t^6 + 1542648801995330874184236*t^5 - 347819053424928336793068*t^4 + 61683475328903338239178*t^3 - 8117056250720937228985*t^2 + 708318740340941449799*t - 30857360406231018655);
C11:=(4*t-5)*(t^2+3*t-6)^6*(9*t^2-28*t+23)^3*(t^4-5*t^3+74*t^2-245*t+223)^3*(4*t^4-9*t^3-t^2+21*t-32)^3*(25*t^4-114*t^3+167*t^2-86*t+20)^3;
     
W11:=[ [F!-27*11^4,54*11^5*43], [F!-27*11^3*131,54*11^4*4973], [] ];
     
     
P13:=[
	t^12+231*t^11+269*t^10-3160*t^9+6022*t^8-9616*t^7+21880*t^6-34102*t^5+28297*t^4-12455*t^3+2876*t^2-243*t+1,
	t^12-9*t^11+29*t^10-40*t^9+22*t^8-16*t^7+40*t^6-22*t^5-23*t^4+25*t^3-4*t^2-3*t+1,
	(t^4-t^3+2*t^2-9*t+3)*(3*t^4-3*t^3-7*t^2+12*t-4)*(4*t^4-4*t^3-5*t^2+3*t-1),
	t^8+235*t^7+1207*t^6+955*t^5+3840*t^4-955*t^3+1207*t^2-235*t+1,
	t^8-5*t^7+7*t^6-5*t^5+5*t^3+7*t^2+5*t+1,
	t^4+7*t^3+20*t^2+19*t+1
	];
Q13:=[
	t^12-512*t^11-13079*t^10-32300*t^9-104792*t^8-111870*t^7-419368*t^6+111870*t^5-104792*t^4+32300*t^3-13079*t^2+512*t+1,
	t^12-8*t^11+25*t^10-44*t^9+40*t^8+18*t^7-40*t^6-18*t^5+40*t^4+44*t^3+25*t^2+8*t+1
	];

Groups["13B.5.2"] := rec<grp| level:=13, zywina:=1, group:=sub<GL(2,Integers(13)) | {[2,0,0,1],[1,1,0,1],[1,0,0,8]}>, jmap := (t^2-t+1)^3*P13[1]^3/((t-1)*t*(t^3-4*t^2+t+1)^13)>;
Groups["13B.5.1"] := rec<grp| level:=13, zywina:=2, group:=sub<GL(2,Integers(13)) | {[8,0,0,1],[1,1,0,1],[1,0,0,2]}>, jmap := (t^2-t+1)^3*P13[2]^3/((t-1)^13*t^13*(t^3-4*t^2+t+1))>;
Groups["13B.5.4"] := rec<grp| level:=13, zywina:=3, group:=sub<GL(2,Integers(13)) | {[2,0,0,2],[1,1,0,1],[1,0,0,8]}>, jmap := -13^4*(t^2-t+1)^3*P13[3]^3/((t^3-4*t^2+t+1)^13*(5*t^3-7*t^2-8*t+5))>;
Groups["13B.4.2"] := rec<grp| level:=13, zywina:=4, group:=sub<GL(2,Integers(13)) | {[2,0,0,1],[1,1,0,1],[1,0,0,4]}>, jmap := (t^4-t^3+5*t^2+t+1)*P13[4]^3/(t*(t^2-3*t-1)^13)>;
Groups["13B.4.1"] := rec<grp| level:=13, zywina:=5, group:=sub<GL(2,Integers(13)) | {[4,0,0,1],[1,1,0,1],[1,0,0,2]}>, jmap := (t^4-t^3+5*t^2+t+1)*P13[5]^3/(t^13*(t^2-3*t-1))>;
Groups["13B"] := rec<grp| level:=13, zywina:=6, group:=sub<GL(2,Integers(13)) | {[2,0,0,1],[1,1,0,1],[1,0,0,2]}>, jmap := (t^2+5*t+13)*P13[6]^3/t>;

// There are three know j-invariants such that G_E(13)=13S4. Then it is not possible to create Groups["13S4"] in the REC format Groups.
//Groups["13S4"] := rec<grp|level:=13,group:=sub<GL(2,Integers(13)) | {[2,0,0,2], [2,0,0,3], [0,-1,1,0], [1,1,-1,1]}>, jmap:=[ (2^4*5*13^4*17^3)/3^13, -(2^12*5^3*11*13^4)/3^13, (2^18*3^3*13^4*127^3*139^3*157^3*283^3*929)/(5^13*61^13) ]>;

W13:=[
	[],[],[],
    	[-27*(t^4-t^3+5*t^2+t+1)^3*P13[4],54*(t^2+1)*(t^4-t^3+5*t^2+t+1)^4*Q13[1]],
	[-27*(t^4-t^3+5*t^2+t+1)^3*P13[5],54*(t^2+1)*(t^4-t^3+5*t^2+t+1)^4*Q13[2]],
	[],[]
	];




Groups["17B.4.2"] := rec<grp|level:=17, zywina:=1, group:=sub<GL(2,Integers(17)) | {[2,0,0,11], [4,0,0,-4], [1,1,0,1]}>,jmap:=F!-17*373^3/2^17>;
Groups["17B.4.6"] := rec<grp|level:=17, zywina:=2, group:=sub<GL(2,Integers(17)) | {[11,0,0,2], [-4,0,0,4], [1,1,0,1]}>, jmap:=F!-17^2*101^3/2>;

Groups["37B.8.1"] := rec<grp|level:=37, zywina:=1, group:=sub<GL(2,Integers(37)) | {[2^3,0,0,1], [1,0,0,2], [1,1,0,1]}>, jmap:=F!-7*11^3>;
Groups["37B.8.2"] := rec<grp|level:=37, zywina:=2, group:=sub<GL(2,Integers(37)) | {[2,0,0,1], [1,0,0,2^3], [1,1,0,1]}>, jmap:=F!-7*137^3*2083^3>;


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



function aux_Q(jE,G)
	case G:
		when "11Nn":
			r:=Roots(Numerator(A11*jE^2+B11*jE+C11));
		 when "11B.10.4":
		 	jG:=Groups[G]`jmap;
		 	if jE eq jG then r:=[[0,1]]; else r:=[]; end if;
		 when "11B.10.5":
		 	jG:=Groups[G]`jmap;
		 	if jE eq jG then r:=[[0,1]]; else r:=[]; end if;
		 when "7Ns.3.1":
		 	jG:=Groups[G]`jmap;
		 	if jE eq jG then r:=[[0,1]]; else r:=[]; end if;
	 	 else
	  		jG:=Groups[G]`jmap;
			r:=Roots(Numerator(jG)-jE*Denominator(jG));
	end case;
	if #r ne 0 then 
		return true, r[1,1];
	  else
		return false, 0;
	end if;
end function;
	
function aux_W(E,p,zyw,i,t0)
	case p:
		when 3: W:=W3[zyw]; d:=-3;
		when 5: W:=W5[zyw]; d:=5;
		when 7: W:=W7[zyw]; d:=-7;
		when 11: W:=W11[zyw]; d:=-11;
		when 13: W:=W13[zyw]; d:=13;
	end case;
	A:=Evaluate(W[1],t0);
	B:=Evaluate(W[2],t0);
	E0:=EllipticCurve([A,B]);
	case i:
		when 1: boo:=IsIsomorphic(E,E0); 
		when 2:	boo:=IsIsomorphic(E,QuadraticTwist(E0,d));
	end case;
	return boo;
end function;


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

function GaloisLabel(E,p)
// Input: E an elliptic curve without CM over Q and p a prime
// Output: label of G_E(p):=Im(rho_{E,p)(Gal(Qbar/Q)) <= GL_2(Fp)
	assert not HasComplexMultiplication(E);

	jE:=jInvariant(E);
	case p:
// p=2
		when 2:
			G3:="2Cn";
			boo3,t3:=aux_Q(jE,G3);
			if boo3 then 
				G1:="2Cs"; 
				boo1,t1:=aux_Q(jE,G1);
				if boo1 then 
					return "2Cs";
				end if;		
				return "2Cn";
			end if;			
			G2:="2B";
			boo2,t2:=aux_Q(jE,G2);
			if boo2 then 
				G1:="2Cs"; 
				boo1,t1:=aux_Q(jE,G1);
				if boo1 then 
					return "2Cs";
				end if;		
				return "2B";
			end if;
			return "GL2(2)";			
// p=3
		when 3:			
			G3:="3B";
			boo3,t3:=aux_Q(jE,G3);
			if boo3 then 
				G1:="3Cs"; 
				boo1,t1:=aux_Q(jE,G1);
				if boo1 then 
					boo11:=aux_W(E,3,1,1,t1);
					boo12:=aux_W(E,3,1,2,t1);
					if boo11 or boo12 then return "3Cs.1.1"; end if;
				  	return "3Cs";
				end if;
				boo31:=aux_W(E,3,3,1,t3);	
				if boo31 then return "3B.1.1"; end if;
				boo32:=aux_W(E,3,3,2,t3);	
				if boo32 then return "3B.1.2"; end if;
				return "3B";
			end if;
			G4:="3Nn";
			boo4,t4:=aux_Q(jE,G4);
			if boo4 then
				G2:="3Ns";
				boo2,t2:=aux_Q(jE,G2);
				if boo2 then 
					G1:="3Cs"; 
				 	boo1,t1:=aux_Q(jE,G1);
				 	if boo1 then 
				 		boo11:=aux_W(E,3,1,1,t1);
						if boo11 then return "3Cs.1.1"; end if;
				  		return "3Cs";
					end if;
					return "3Ns";
				end if;
				return "3Nn"; 	
			end if;
			return "GL2(3)";
// p=5
		when 5:
			G8:="5B";
			boo8,t8:=aux_Q(jE,G8);
			if boo8 then 
				G2:="5Cs"; 
				boo2,t2:=aux_Q(jE,G2);
				if boo2 then
					G1:="5Cs.4.1";
					boo1,t1:=aux_Q(jE,G1); 
					if boo1 then 
						boo11:=aux_W(E,5,1,1,t1);
						if boo11 then return "5Cs.1.1"; end if;
						boo12:=aux_W(E,5,1,2,t1);
						if boo12 then return "5Cs.1.3"; end if;					
						return "5Cs.4.1";
					end if;
					return "5Cs";
				end if;
				G6:="5B.4.1"; 
				boo6,t6:=aux_Q(jE,G6);
				if boo6 then
					G1:="5Cs.4.1";
					boo1,t1:=aux_Q(jE,G1); 
					if boo1 then 
						boo11:=aux_W(E,5,1,1,t1);
						if boo11 then return "5Cs.1.1"; end if;
						boo12:=aux_W(E,5,1,2,t1);
						if boo12 then return "5Cs.1.3"; end if;					
						return "5Cs.4.1";
					end if;
					boo61:=aux_W(E,5,6,1,t6);
					if boo61 then return "5B.1.1"; end if;
					boo62:=aux_W(E,5,6,2,t6);
					if boo62 then return "5B.1.4"; end if;	
					return "5B.4.1";
				end if;
				G5:="5B.4.2"; 
				boo5,t5:=aux_Q(jE,G5);
				if boo5 then
					G1:="5Cs.4.1";
					boo1,t1:=aux_Q(jE,G1); 
					if boo1 then 
						boo11:=aux_W(E,5,1,1,t1);
						if boo11 then return "5Cs.1.1"; end if;
						boo12:=aux_W(E,5,1,2,t1);
						if boo12 then return "5Cs.1.3"; end if;					
						return "5Cs.4.1";
					end if;
					boo51:=aux_W(E,5,5,1,t5);
					if boo51 then return "5B.1.2"; end if;
					boo52:=aux_W(E,5,5,2,t5);
					if boo52 then return "5B.1.3"; end if;
					return "5B.4.2";	
				end if;
				return "5B";
			end if;
			G9:="5S4";
			boo9,t9:=aux_Q(jE,G9);
			if boo9 then 
				G4:="5Ns";
				boo4,t4:=aux_Q(jE,G4);
				if boo4 then
					G2:="5Cs"; 
					boo2,t2:=aux_Q(jE,G2);
					if boo2 then
						G1:="5Cs.4.1";
						boo1,t1:=aux_Q(jE,G1); 
						if boo1 then 	
							boo11:=aux_W(E,5,1,1,t1);
							if boo11 then return "5Cs.1.1"; end if;
							boo12:=aux_W(E,5,1,2,t1);
							if boo12 then return "5Cs.1.3"; end if;					
							return "5Cs.4.1";
						end if;		
						return "5Cs";
					end if;
					G3:="5Ns.2.1";
					boo3,t3:=aux_Q(jE,G3);
					if boo3 then return "5Ns.2.1"; end if;
					return "5Ns";
				end if;
				return "5S4";
			end if;
			G7:="5Nn";
			boo7,t7:=aux_Q(jE,G7);
			if boo7 then
				G3:="5Ns.2.1";
				boo3,t3:=aux_Q(jE,G3);
				if boo3 then return "5Ns.2.1"; end if;
				return "5Nn";
			end if;	
			return "GL2(5)"; 
// p=7																					
		when 7:			
			G7:="7B";
			boo7,t7:=aux_Q(jE,G7);
			if boo7 then 
				G3:="7B.6.1"; 
				boo3,t3:=aux_Q(jE,G3);
				if boo3 then 
					boo31:=aux_W(E,7,3,1,t3);
					if boo31 then return "7B.1.1"; end if;
					boo32:=aux_W(E,7,3,2,t3);
					if boo32 then return "7B.1.6"; end if;	
					return "7B.6.1";				
				end if;
				G4:="7B.6.3"; 
				boo4,t4:=aux_Q(jE,G4);				
				if boo4 then 
					boo41:=aux_W(E,7,4,1,t4);
					if boo41 then return "7B.1.3"; end if;
					boo42:=aux_W(E,7,4,2,t4);
					if boo42 then return "7B.1.4"; end if;	
					return "7B.6.3";				
				end if;				
				G5:="7B.6.2"; 	
				boo5,t5:=aux_Q(jE,G5);			
				if boo5 then 
					boo51:=aux_W(E,7,5,1,t5);
					if boo51 then return "7B.1.5"; end if;
					boo52:=aux_W(E,7,5,2,t5);
					if boo52 then return "7B.1.2"; end if;	
					return "7B.6.2";				
				end if;				
				boo71:=aux_W(E,7,7,1,t7);
				if boo71 then return "7B.2.3"; end if;
				boo72:=aux_W(E,7,7,2,t7);
				if boo72 then return "7B.2.1"; end if;		
				return "7B";
			end if;
			G6:="7Nn";
			boo6,t6:=aux_Q(jE,G6);
			if boo6 then return "7Nn"; end if;
			G2:="7Ns";
			boo2,t2:=aux_Q(jE,G2);
			if boo2 then 
				G1:="7Ns.3.1";
				boo1,t1:=aux_Q(jE,G1);
				if boo1 then 
					boo11:=aux_W(E,7,1,1,t1);
					boo12:=aux_W(E,7,1,2,t1);
					if boo11 or boo12 then return "7Ns.2.1"; end if;
					return "7Ns.3.1";
				end if;
				return "7Ns";
			end if;
			return "GL2(7)"; 
// p=11
		when 11:			
			G3:="11Nn";
			boo3,t3:=aux_Q(jE,G3);
			if boo3 then 
				return "11Nn";
			end if;
			G1:="11B.10.4";
			boo1,t1:=aux_Q(jE,G1);
			if boo1 then					
				boo11:=aux_W(E,11,1,1,t1);	
				if boo11 then return "11B.1.4"; end if;
				boo12:=aux_W(E,11,1,2,t1);	
				if boo12 then return "11B.1.7"; end if;
				return "11B.10.4";
			end if;
			G2:="11B.10.5";
			boo2,t2:=aux_Q(jE,G2);
			if boo2 then				
				boo21:=aux_W(E,11,2,1,t1);	
				if boo21 then return "11B.1.5"; end if;
				boo22:=aux_W(E,11,2,2,t1);	
				if boo22 then return "11B.1.6"; end if;
				return "11B.10.5";
			end if;
			return "GL2(11)";
// p=13
		when 13:
			G6:="13B";
			boo6,t6:=aux_Q(jE,G6);
			if boo6 then 
				G1:="13B.5.2";
				boo1,t1:=aux_Q(jE,G1);
				if boo1	then return "13B.5.2"; end if;
				G2:="13B.5.1";
				boo2,t2:=aux_Q(jE,G2);
				if boo2	then return "13B.5.1"; end if;							
				G3:="13B.5.4";
				boo3,t3:=aux_Q(jE,G3);
				if boo3	then return "13B.5.4"; end if;			
				G4:="13B.4.2";
				boo4,t4:=aux_Q(jE,G4);
				if boo4	then					
					boo41:=aux_W(E,13,4,1,t4);	
					if boo41 then return "13B.3.2"; end if;
					boo42:=aux_W(E,13,4,2,t4);	
					if boo42 then return "13B.3.7"; end if;
					return "13B.4.2";
				end if;
				G5:="13B.4.1";
				boo5,t5:=aux_Q(jE,G5);
				if boo5	then					
					boo51:=aux_W(E,13,5,1,t5);	
					if boo51 then return "13B.3.1"; end if;
					boo52:=aux_W(E,13,5,2,t5);	
					if boo52 then return "13B.3.4"; end if;
						return "13B.4.1";
					end if;					
				return "13B";
			end if;
			if jE in {(2^4*5*13^4*17^3)/3^13, -(2^12*5^3*11*13^4)/3^13, (2^18*3^3*13^4*127^3*139^3*157^3*283^3*929)/(5^13*61^13)} then return "13S4"; end if;					
			// Denominator of jE has 13th-power (Banwait-Cremona: Tetraedral...
			// Dogra (STNB19): Determine X_S4(13)(Q)
			if aux_IsSurjective(E,13) then return "GL2(13)"; end if;
			// Zywina: section 6.3: to determine GL2(13)
			// Sutherland: Compute explicitly G_E(13) 						
			return "GL2(13)";
// p=17
		when 17:
			case jE: 
				when -17*373^3/2^17: return "17B.4.2";
				when -17^2*101^3/2: return "17B.4.6";
			  else
			if aux_IsSurjective(E,17) then return "GL2(17)"; end if;
			end case;
// p=37
		when 37:
			case jE: 
				when -7*11^3: return "37B.8.1";
				when -7*137^3*2083^3: return "37B.8.2";
			  else
			if aux_IsSurjective(E,37) then return "GL2(37)"; end if;
			end case;						  					
	end case;
	if aux_IsSurjective(E,p) then return "GL2(" cat IntegerToString(p) cat ")"; end if;
	return "Counterexample to Uniformity Conjecture";
end function;





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

function ExceptionalGaloisLabels(E)
// Input: E an elliptic curve without CM over Q 
// Output: List of labels of G_E(p):=Im(rho_{E,p)(Gal(Qbar/Q)) < GL_2(Fp) for p exceptional prime
	assert not HasComplexMultiplication(E);
	SE:=ExceptionalSet(E); 
	L:=[];
	for p in SE do
		GEp:=GaloisLabel(E,p);
		if GEp ne "GL2(" cat IntegerToString(p) cat ")" then 
			L:=Append(L,GEp);
		end if; 
	end for;
	if booprintGal then printf "Exceptional Galois Representation (non-CM):%o\n",L; end if;
	return L;
end function;



////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// Enrique Gonzalez-Jimenez & Filip Najman,"Growth of torsion groups of elliptic curves upon base changes". arXiv: 1609.02515
//
// Section 6
//
// http://verso.mat.uam.es/~enrique.gonzalez.jimenez/research/tables/growth/RQ.txt
// 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

function C(p)
// Input: prime p
// Output: the set A with the property that it is the smallest set values a such that p is in R_Q(d) if and only if d is divisible by some a in A.

CM:={3,4,7,8,11,19,43,67,163};
  if p in {2,3,5,7} then return [1]; end if;
  if p in {11} then return [5]; end if;
  if p in {13} then return [3,4]; end if;
  if p in {17} then return [8]; end if;
  if p in {37} then return [12]; end if;
  if p in {19,43,67,163} then return [(p-1)/2]; end if;
  if p mod 3 eq 1 then return [2*(p-1)]; end if;
  if p mod 3 eq 2 and not &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [2*(p-1)]; end if;
  if p mod 9 in {2,5} and &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [(p^2-1)/3]; end if;
  if p mod 9 eq 8 and &and[LegendreSymbol(-D,p) eq -1 : D in CM] then 
    return [p^2-1,(p^2-1)/3]; 
  end if;
end function;

function RQ(d)
// Given a positive integer d, 
// the function returns the set of primes p for which 
// there exists a number field K with [K:Q]=d 
// and an elliptic curve E/Q , 
// such that p divides the order of the torsion subgroup of E(K)

  bound:=Integers()!Max((163-1)/2,2*d+1);
  R:={2,3,5,7};
  if d le 2 then return R; end if;
  for p in PrimesInInterval(2,bound) do Cp:=C(p); 
    for cp in Cp do
      if IsDivisibleBy(d,Integers()!cp) then R:=R join {p}; end if;
     end for; 
  end for;
  return R;

end function;


////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// Enrique Gonzalez-Jimenez & Filip Najman,"Growth of torsion groups of elliptic curves upon base changes". arXiv: 1609.02515
//
// Tables 1 & 2
//
// http://verso.mat.uam.es/~enrique.gonzalez.jimenez/research/tables/growth/RQ.txt
// 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// non-CM


function ExcepcionalGalrepDegrees(galrep)
// Given the label(=galrep) of a non-surjective mod p galois rep 
// attached to an elliptic curve E over Q without CM 
// the function returns {[Q(P):Q] , ord(P)=p ]}, [Q(E[p]):Q]

  case galrep:
    when "2Cs": return { 1 },1;
    when "2B": return { 1, 2 },2;
    when "2Cn": return { 3 },3;
    when "3Cs.1.1": return { 1, 2 },2;
    when "3Cs": return { 2, 4 },4;
    when "3B.1.1": return { 1, 6 },6;
    when "3B.1.2": return { 2, 3 },6;
    when "3Ns": return { 4 },8;
    when "3B": return { 2, 6 },12;
    when "3Nn": return { 8 },16;
    when "5Cs.1.1": return { 1, 4 },4;
    when "5Cs.1.3": return { 2, 4 },4;
    when "5Cs.4.1": return { 2, 4, 8 },8;
    when "5Ns.2.1": return { 8, 16 },16;
    when "5Cs": return { 4, 16 },16;
    when "5B.1.1": return { 1, 20 },20;
    when "5B.1.2": return { 4, 5 },20;
    when "5B.1.3": return { 4, 10 },20;
    when "5B.1.4": return { 2, 20 },20;
    when "5Ns": return { 8, 16 },32;
    when "5B.4.1": return { 2, 20 },40;
    when "5B.4.2": return { 4, 10 },40;
    when "5Nn": return { 24 },48;
    when "5B": return { 4, 20 },80;
    when "5S4": return { 24 },96;
    when "7Ns.2.1": return { 6, 9, 18 },18;
    when "7Ns.3.1": return { 12, 18 },36;
    when "7B.1.1": return { 1, 42 },42;
    when "7B.1.2": return { 3, 42 },42;
    when "7B.1.3": return { 6, 7 },42;
    when "7B.1.4": return { 3, 14 },42;
    when "7B.1.5": return { 6, 21 },42;
    when "7B.1.6": return { 2, 21 },42;
    when "7Ns": return { 12, 36 },72;
    when "7B.6.1": return { 2, 42 },84;
    when "7B.6.2": return { 6, 42 },84;
    when "7B.6.3": return { 6, 14 },84;
    when "7Nn": return { 48 },96;
    when "7B.2.1": return { 3, 42 },126;
    when "7B.2.3": return { 6, 21 },126;
    when "7B": return { 6, 42 },252;
    when "11B.1.4": return { 5, 110 },110;
    when "11B.1.5": return { 5, 110 },110;
    when "11B.1.6": return { 10, 55 },110;
    when "11B.1.7": return { 10, 55 },110;
    when "11B.10.4": return { 10, 110 },220;
    when "11B.10.5": return { 10, 110 },220;
    when "11Nn": return { 120 },240;
    when "13S4": return { 72, 96 },288;
    when "13B.3.1": return { 3, 156 },468;
    when "13B.3.2": return { 12, 39 },468;
    when "13B.3.4": return { 6, 156 },468;
    when "13B.3.7": return { 12, 78 },468;
    when "13B.5.1": return { 4, 156 },624;
    when "13B.5.2": return { 12, 52 },624;
    when "13B.5.4": return { 12, 156 },624;
    when "13B.4.1": return { 6, 156 },936;
    when "13B.4.2": return { 12, 78 },936;
    when "13B": return { 12, 156 },1872;
    when "17B.4.2": return { 8, 272 },1088;
    when "17B.4.6": return { 16, 136 },1088;
    when "37B.8.1": return { 12, 1332 },15984;
    when "37B.8.2": return { 36, 444 },15984;
    when "37GL": return { 1368 },1822176;
  end case;
end function;

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

function aux_number(label)
// Given a label(=string) of the form, for example 23ab2, the function return the number before the first letter, in this case 23.  
// It will be used in:
// prime of the label for the galrep
// level of the Cremona label of an elliptic curve over Q.
  S:=Eltseq(label);
  P:=[];
  c:=1;
  while StringToCode(S[c]) in [48..57] do // check that the string is a number 0-9
    P:=Append(P,S[c]); 
    c:=c+1;
  end while; 
  return StringToInteger(&cat[s : s in P]);
end function;



function PrimesNoCM(E,d)
// Input: E an elliptic curve over Q WITHOUT CM and a positive integer d
// Output: Set of pairs <p,boop1,boop2> such that p is a primes such that 
//            there exists a point P of order p such that [Q(P):Q] divides d 
//            and
//            	 boop1 is true if E(Q)[p] \ne {0}, otherwise false
//               boop2 is true if [Q(E[p]) : Q] divides d, otherwise false

//  label,galreps:=aux_split(S);
//  E:=EllipticCurve(label);
	galreps:=ExceptionalGaloisLabels(E);
	//printf "Non-CM: Non-surjective images:%o\n", galreps;
	assert not HasComplexMultiplication(E);
	A:=<<aux_number(galrep),galrep> : galrep in galreps>;
	P1:={p[1] : p in A}; 
//	P2:=RQpre(d) diff P1;
	P2:=RQ(d) diff P1;
	L:=[];
	for galrep in A do
    	l:=galrep[1];
    	lab:=galrep[2];
    	dP,dE:=ExcepcionalGalrepDegrees(lab);
		booQ:=1 in dP;
		booFull:=IsDivisibleBy(d,dE);
    	bool:=&or[IsDivisibleBy(d,dl) : dl in dP];
		if bool then 		
			L:=Append(L,<l,booQ,booFull>);
		end if;
	end for;
	for l in P2 do
		dl:=l^2-1;
    	dE:=l*(l-1)*(l^2-1);
    	bool:=IsDivisibleBy(d,dl);
    	booQ:=false; // No points of order l defined over Q. 
    	booFull:=IsDivisibleBy(d,dE); 
 		if bool then 		
			L:=Append(L,<l,booQ,booFull>);
		end if;
	end for;
	return Sort(L);
end function;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
//  D. Zywina, "On the possible images of the mod l representations associated to elliptic curves over Q". arXiv:1508.07660
//
// Section 1.8
//
// Enrique Gonzalez-Jimenez & Filip Najman,"Growth of torsion groups of elliptic curves upon base changes". arXiv: 1609.02515
//
// 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// CM


function CnoCM(p)
// Given a prime p, the function returns the set A with the property that it is the smallest set values a such that p is in R^CM_Q(d) if and only if d is divisible by some a in A.

CM:={3,4,7,8,11,19,43,67,163};
  if p in {2,3} then return [1]; end if;
  if p in {7} then return [3]; end if;
  if p in {11} then return [5]; end if;
  if p in {5} then return [8]; end if;
  if p in {13} then return [24]; end if;
  if p in {19,43,67,163} then return [(p-1)/2]; end if;
  if p mod 3 eq 1 then return [2*(p-1)]; end if;
  if p mod 3 eq 2 and not &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [2*(p-1)]; end if;
  if p mod 9 in {2,5} and &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [(p^2-1)/3]; end if;
  if p mod 9 eq 8 and &and[LegendreSymbol(-D,p) eq -1 : D in CM] then return [p^2-1,(p^2-1)/3]; end if;
end function;

function RCMQ(d)
// Given a positive integer d, 
// the function returns the set of primes p for which 
// there exists a number field K with [K:Q]=d 
// and an elliptic curve E/Q with CM, 
// such that p divides the order of the torsion subgroup of E(K)

  bound:=Integers()!Max((163-1)/2,2*d+1);
  R:={2,3};
  if d le 2 then return R; end if;
  for p in PrimesInInterval(2,bound) do Cp:=CnoCM(p); 
    for cp in Cp do
      if IsDivisibleBy(d,Integers()!cp) then R:=R join {p}; end if;
     end for; 
  end for;
  return R;

end function;

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


function DegCM_2(D)
  if D in {-12,-16,-7,-28,-8} then dv:=[1,2]; end if;
  if D in {-43,-27,-11,-163,-19,-67} then dv:=[3]; end if;
  if D eq -4 then dv:=[1,2]; end if;
  if D eq -3 then dv:=[1,2,3]; end if;
  return dv;
end function;

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

function DegsCM_jne0(l,D)
// l odd
// Prop 1.14 Zywina:
// j ne 0
  assert IsPrime(l) and D ne -3 and IsOdd(l);

  case IsDivisibleBy(D,l):
    when false: 
      case IsSquare(FiniteField(l)!D):
        when true: dv:=[2*(l-1),(l-1)^2];
        when false: dv:=[l^2-1];
     end case;
   when true:  
     dv:=[(l-1)/2,l-1,(l-1)*l/2,(l-1)*l];
  end case;
  return dv;
end function;

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

function DegsCM_jeq0(l)
// l odd
// Prop 1.16 Zywina:
// j=0
  assert IsPrime(l) and IsOdd(l);

  case l mod 9:
    when 1: dv:=[2*(l-1),(l-1)^2];
    when 8: dv:=[l^2-1];
    when 4: dv:=[2*(l-1),(l-1)^2/3,2*(l-1)^2/3,(l-1)^2];           
    when 7: dv:=[2*(l-1),(l-1)^2/3,2*(l-1)^2/3,(l-1)^2];
    when 2: dv:=[(l^2-1)/3,2*(l^2-1)/3,l^2-1];
    when 5: dv:=[(l^2-1)/3,2*(l^2-1)/3,l^2-1];                         
    when 3: dv:=[1,2,3,4,6];
   end case;
  return dv;
end function;

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

function DegsCM(l,D)
 if l eq 2 then return {Integers()!d: d in DegCM_2(D)}; end if;
 if D eq -3 then
    return {Integers()!d: d in DegsCM_jeq0(l)};
   else
    return {Integers()!d: d in DegsCM_jne0(l,D)};
  end if;
end function;

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

function DegsCMall(l)
  return &join[DegsCM(l,D) : D in [-3,-12,-27,-4,-16,-7,-28,-8,-11,-19,-43,-67,-163]];
end function;

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

function CM_EllipticCurve(D,d)
  case D: 
    when -3:   s:=[ 0, 16 ]; E:=EllipticCurve([0,16*d]); 
    when -12:  s:=[ -15, 22 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -27:  s:=[ -480, 4048 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -4:   s:=[ 1, 0 ]; E:=EllipticCurve([d,0]); 
    when -16:  s:=[ -11, 14 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -7:   s:=[ -1715, 33614 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -28:  s:=[ -29155, 1915998 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -8:   s:=[ -4320, 96768 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -11:  s:=[ -9504, 365904 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -19:  s:=[ -608, 5776 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -43:  s:=[ -13760, 621264 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -67:  s:=[ -117920, 15585808 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -163: s:=[ -34790720, 78984748304 ]; E:=QuadraticTwist(EllipticCurve(s),d);
  end case;
  
  return E;
end function;

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

function IsCube(d)
  Num:=Numerator(d); 
  Den:=Denominator(d);
  if &and[f[2] mod 3 eq 0 : f in Factorization(Num)] and &and[f[2] mod 3 eq 0 : f in Factorization(Den)] then
    return true;
   else
    return false;
  end if;
end function;

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

function DegCM_2D(E)
  boo,D:=HasComplexMultiplication(E);
  assert boo;
  if D in {-12,-16,-7,-28,-8} then dv:=[1,2]; dE:=2; end if;
  if D in {-43,-27,-11,-163,-19,-67} then dv:=[3]; dE:=6; end if;
  if D eq -4 then 
    c:=cInvariants(E);
    c4:=c[1];c6:=c[2];
    assert c6 eq 0;
    d:=27*c4;
    if IsSquare(d) then dv:=[1]; dE:=1; else dv:=[1,2]; dE:=2; end if;
  end if;
  if D eq -3 then 
    c:=cInvariants(E);
    c4:=c[1];c6:=c[2];
    assert c4 eq 0;
    d:=-54*c6;
    if IsCube(d) then
      dv:=[1,2]; dE:=2;
     else
      dv:=[3]; dE:=6;
    end if;
  end if;
  return dv,dE;
end function;

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

function GECM_2j(E)
  boo,D:=HasComplexMultiplication(E);
  assert boo;
  jE:=jInvariant(E);
  if jE in {2^4*3^3*5^3,2^3*3^3*11^3,-3^3*5^3,3^3*5^3*17^3,2^6*5^3} then dv:=[1,2]; dE:=2; end if;
  if jE in {-2^(15)*3*5^3,-2^(15),-2^(15)*3^3,-2^(18)*3^3*5^3,-2^(15)*3^3*5^3*11^3,-2^(18)*3^3*5^3*23^3*29^3} then dv:=[3]; dE:=6; end if;
  if jE eq 1728 then 
    c:=cInvariants(E);
    c4:=c[1];c6:=c[2];
    assert c6 eq 0;
    d:=27*c4;
    if IsSquare(d) then dv:=[1]; dE:=1; else dv:=[1,2]; dE:=2; end if;
  end if;
  if jE eq 0 then 
    c:=cInvariants(E);
    c4:=c[1];c6:=c[2];
    assert c4 eq 0;
    d:=-54*c6;
    if IsCube(d) then
      dv:=[1,2]; dE:=2;
     else
      dv:=[3]; dE:=6;
    end if;
  end if;
  return dv,dE;
end function;

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

function DegCM_jne0_lodd(E,l)
  boo,D:=HasComplexMultiplication(E);
  assert boo and (jInvariant(E) ne 0) and IsPrime(l) and IsOdd(l);
    
  Fl:=FiniteField(l);
  Fl0:=[a : a in Fl | not IsZero(a)];
  Fl2:=[a : a in Fl0 | IsSquare(a)];
  if IsDivisibleBy(D,l) then 
     if IsIsomorphic(E,CM_EllipticCurve(D,1)) then   
        dv:=[(l-1)/2,(l-1)*l]; dE:=l*(l-1);
        return dv,dE;
       else
           if IsIsomorphic(E,CM_EllipticCurve(D,-l)) then   
              dv:=[(l-1),(l-1)*l/2]; dE:=l*(l-1);
              return dv, dE;
            else 
               dv:=[(l-1),(l-1)*l]; dE:=2*l*(l-1);
               return dv, dE;
           end if;
      end if;
    else
     if IsSquare(Fl!D) then 
        dv:=[2*(l-1),(l-1)^2]; dE:=2*(l-1)^2;
        return dv, dE;
      else
        dv:=[l^2-1]; dE:=2*(l^2-1);
        return dv, dE;
     end if;
   end if;
end function;

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

function DegCM_j0_lodd(E,l)
  boo,D:=HasComplexMultiplication(E);
  assert boo and (jInvariant(E) eq 0) and (D eq -3) and IsPrime(l) and IsOdd(l);
    
  case l mod 9:
    when 1:
           dv:=[2*(l-1),(l-1)^2]; dE:=2*(l-1)^2;
           return dv, dE;
    when 8:
           dv:=[(l^2-1)];  dE:=2*(l^2-1);
           return dv,dE;
    when 4:
           ee:=Integers()!((l-1)/3);
           e:=ee mod 3; 
           assert e in {1,2};
           EE:=EllipticCurve([0,16*l^e]);/////////////////  
           if IsQuadraticTwist(E,EE) then
             dv:=[2*(l-1),(l-1)^2/3,2*(l-1)^2/3]; dE:=2/3*(l-1)^2;
             return dv, dE;
            else
             dv:=[2*(l-1),(l-1)^2]; dE:=2*(l-1)^2;
             return dv, dE;
           end if;
    when 7:
           ee:=Integers()!((l-1)/3);
           e:=ee mod 3;  
           assert e in {1,2};
           EE:=EllipticCurve([0,16*l^e]);/////////////////  
           if IsQuadraticTwist(E,EE) then
             dv:=[2*(l-1),(l-1)^2/3,2*(l-1)^2/3];  dE:=2/3*(l-1)^2;
             return dv, dE;
            else
             dv:=[2*(l-1),(l-1)^2]; dE:=2*(l-1)^2;
             return dv, dE;
           end if;
    when 2:
           ee:=Integers()!(-(l+1)/3);
           e:=ee mod 3; 
           assert e in {1,2};
           EE:=EllipticCurve([0,16*l^e]);/////////////////  
           if IsQuadraticTwist(E,EE) then
             dv:=[(l^2-1)/3,2*(l^2-1)/3]; dE:=2/3*(l^2-1);
             return dv, dE;
            else
             dv:=[(l^2-1)]; dE:=2*(l^2-1);
             return dv, dE;
           end if;                             
    when 5:
           ee:=Integers()!(-(l+1)/3);
           e:=ee mod 3; 
           assert e in {1,2};
           EE:=EllipticCurve([0,16*l^e]);/////////////////  
           if IsQuadraticTwist(E,EE) then
             dv:=[(l^2-1)/3,2*(l^2-1)/3];dE:=2/3*(l^2-1);
             return dv, dE;
            else
             dv:=[(l^2-1)];dE:=2*(l^2-1);
             return dv, dE;
           end if;                             
    when 3:
           c:=cInvariants(E);
           c4:=c[1];c6:=c[2];
           assert c4 eq 0;
           d:=-54*c6;

           if &or[IsSquare(d), IsSquare(-3*d)] and IsCube(-4*d) then 
             dv:=[1,2]; dE:=2;
             return dv, dE;
           end if;

          if not IsSquare(d) and not IsSquare(-3*d) and IsCube(-4*d) then 
             dv:=[2,4]; dE:=4;
             return dv, dE;
           end if;

          if IsSquare(d) and not IsCube(-4*d) then 
             dv:=[1,6]; dE:=6;
             return dv, dE;
           end if;

          if IsSquare(-3*d) and not IsCube(-4*d) then 
             dv:=[2,3]; dE:=6;
             return dv, dE;
           end if;

          if not IsSquare(d) and not IsSquare(-3*d) and not IsCube(-4*d) then 
             dv:=[2,6]; dE:=12;
             return dv, dE;
           end if;
  
    end case;
end function;

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

function DegreesCM(E,l)
  boo,D:=HasComplexMultiplication(E);
  assert boo;
  if l eq 2 then dv,dE:=GECM_2j(E); return {Integers()!d: d in dv}, Integers()!dE; end if;
  if jInvariant(E) eq 0 then 
    dv,dE:=DegCM_j0_lodd(E,l);
   else
    dv,dE:=DegCM_jne0_lodd(E,l);
  end if;
  return {Integers()!d: d in dv}, Integers()!dE; 
end function;


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


function PrimesCM(E,d)
// Input: E an elliptic curve over Q WITH CM and a positive integer d
// Output: Set of pairs <p,boop1,boop2> such that p is a primes such that 
//            there exists a point P of order p such that [Q(P):Q] divides d 
//            and
//            	 boop1 is true if E(Q)[p] \ne {0}, otherwise false
//               boop2 is true if [Q(E[p]) : Q] divides d, otherwise false

//  label,galreps:=aux_split(S);
//  E:=EllipticCurve(label);
	assert HasComplexMultiplication(E);
	L:=[];
	for l in RCMQ(d) do
		dP,dE:=DegreesCM(E,l);
		booQ:=1 in dP;
		booFull:=IsDivisibleBy(d,dE);
    	bool:=&or[IsDivisibleBy(d,dl) : dl in dP];
		if bool then 		
			L:=Append(L,<l,booQ,booFull>);
		end if;
	end for;
	return L;
end function;



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

function R(E,d)
// Input: E an elliptic curve defined over Q and d a positive integer
// Output: Set of primes where there exists a number field K of degree dividing d such that the p-power torsion growth over K.
//         a second output is given in some cases: a list of points and number field (TODO)
	REd:=[];
	if #[p : p in PrimeFactors(d) | p lt 11] eq 0 then return [],[]; end if;
	booCM,D:=HasComplexMultiplication(E);
	if booCM then 
		RE:=PrimesCM(E,d);
	  else
	  	RE:=PrimesNoCM(E,d);
	end if;
	for sp in RE do
		p:=sp[1];
		booQ:=sp[2];
		booFull:=sp[3];
		case p:
			when 2:
				T2:=TwoTorsionSubgroup(E);
				if booFull and #T2 ne 4 then
					REd:=Append(REd,p);
				  else
					if not booQ then 
						REd:=Append(REd,p);
					  else
					  	t2:=AbelianInvariants(T2);
					  	if #t2 eq 1 then
					  		k:=Valuation(t2[1],2);
					  		if GCD(p*(p-1),d) ne 1 and #PolDegree0(E,2,k+1,d) ne 0 then
								REd:=Append(REd,p); 
							end if;
					  	  else		  	
					  	  	k:=Valuation(t2[2],2);
							if GCD(p*(p-1),d) ne 1 then
								if #PolDegree0(E,2,2,d) ne 0 then
									REd:=Append(REd,p); 
								  else
								    if k gt 1 then 
								    	if	#PolDegree0(E,2,k+1,d) ne 0 then
										REd:=Append(REd,p);  
										end if;
									end if;
								end if;
							end if;
						end if;
					end if;
				end if;									 			  
			when 3:
				if booFull then 
					REd:=Append(REd,p); 
				  else
					if not booQ then 
						REd:=Append(REd,p);
					  else
						T:=TorsionSubgroup(E);
						k:=Valuation(#T,3);
						if GCD(p*(p-1),d) ne 1 and #PolDegree0(E,3,k+1,d) ne 0 then
								REd:=Append(REd,p); 
						end if;
					end if;
				end if;						
			when 5:
				if booFull then 
					REd:=Append(REd,p); 
				  else
					if not booQ then 
						REd:=Append(REd,p);
					  else
						if GCD(p*(p-1),d) ne 1 and #PolDegree0(E,5,2,d) ne 0 then
								REd:=Append(REd,p); 
						end if;
					end if;
				end if;				
			when 7:
				if booFull then 
					REd:=Append(REd,p); 
				  else
					if not booQ then 
						REd:=Append(REd,p);
					  else
						if GCD(p*(p-1),d) ne 1 and #PolDegree0(E,7,2,d) ne 0 and d ge 42 then
								REd:=Append(REd,p); 
						end if;
					end if;
				end if;		
			else REd:=Append(REd,p);						
		end case;		
	end for;
	return REd;
end function;




// TEST

function aux_space(S)
// Given a string of characters S the function returns the places where the space=" " are
// For example for each line S in the file from the Cremona database

  L:=[];
  for k in [1..#S] do
    if S[k] eq " " then
      L:=Append(L,k);
    end if; 
  end for;
  L:=Append(L,#S);
  return L;
end function;

function aux_split(A)
// Given a string of characters A corresponding to an elliptic curve E over Q and the galois representation for each excepcional primes the function returns Cremona label and Sutherland Galois representations label 
// For example for each line A (=elliptic curve) in the file  from the Cremona database  
  S:=Eltseq(A);
  spaces:=aux_space(S);
  if #spaces eq 1 then return A,[]; end if;
  label:=Substring(A,1,spaces[1]-1);
  galrep:=[];
  for k in [1..#spaces-2] do 
    star:=spaces[k]+1;
    fin:=spaces[k+1]-spaces[k]-1;
    galrep:=Append(galrep,Substring(A,star,fin));
  end for;
  k:=#spaces-1;
  star:=spaces[k]+1;
  fin:=spaces[k+1];
  galrep:=Append(galrep,Substring(A,star,fin));
  return label,galrep; 
end function;







/*
// TEST 1

	data:="galrep.000000-399999";
	//data:="galrep.00000-09999";
	SS:=Split(Read(data));
	for S in SS do
		label,galreps:=aux_split(S);
		E:=EllipticCurve(label);
		//assert not HasComplexMultiplication(E);
		if not HasComplexMultiplication(E) then 
			A:=<<aux_number(galrep),galrep> : galrep in galreps>;
			P1:={p[1] : p in A}; 
			label,galreps;
			for a in A do 
				p:=a[1];
				GEp:=GaloisLabel(E,p);
				if GEp ne a[2] then print S; break; end if;
			end for;
		end if;
	end for;

// TEST 2

	data:="galrep.000000-399999";
	SS:=Split(Read(data));
	for S in SS do
		label,galreps:=aux_split(S);
		E:=EllipticCurve(label);
		//assert not HasComplexMultiplication(E);
		if not HasComplexMultiplication(E) then 
			Ex:=ExceptionalGaloisLabels(E);
			label,galreps,Ex;
			if Ex ne galreps then break; end if;
		end if;
	end for;
	
// TEST 3

	data:="galrep.000000-399999";
	SS:=Split(Read(data));
	for S in SS do
		label,galreps:=aux_split(S);
		E:=EllipticCurve(label);
		if not HasComplexMultiplication(E) then 
			label,galreps;
			if not &and[aux_IsSurjective(E,p) : p in PrimesInInterval(13,100)] then break; end if;
		end if;
	end for;

*/
