// C2XC30

P< x >:=PolynomialAlgebra(Rationals());

E:=EllipticCurve([1,1,1,-3,1]); //j-Invariant equals -121945/32

f2:=DivisionPolynomial(E,2);
f3:=DivisionPolynomial(E,3);
f5:=DivisionPolynomial(E,5);
f6:=DivisionPolynomial(E,6);
f10:=DivisionPolynomial(E,10);
f15:=DivisionPolynomial(E,15);
f30:=DivisionPolynomial(E,30);

g6:=f6 div f2;
g6:=g6 div f3;

g10:=f10 div f2;
g10:=g10 div f5; 

g15:=f15 div f3;
g15:=g15 div f5;

g30:=f30 div f2;
g30:=g30 div f3;
g30:=g30 div f5;
g30:=g30 div g6;
g30:=g30 div g10;
g30:=g30 div g15; //g30 is now a primitive 30th division polynomial attached to E.

Factorisation(g30: DegreeLimit:=6);

//This outputs polynomials f and g defined below.

f:=x^6-40*x^5+45*x^4-120*x^3+75*x^2-25;
g:=x^6+10*x^5+25*x^4-20*x^3-25*x^2+50*x-25;

Order(GaloisGroup(f));
Order(GaloisGroup(g));
//Splitting fields of f and g  are fields K1,K2 such that [K_1:Q]=[K_2:Q]=12.

E:=EllipticCurve([1,1,1,549,-2202]);
f2:=DivisionPolynomial(E,2);
f3:=DivisionPolynomial(E,3);
f5:=DivisionPolynomial(E,5);
f6:=DivisionPolynomial(E,6);
f10:=DivisionPolynomial(E,10);
f15:=DivisionPolynomial(E,15);
f30:=DivisionPolynomial(E,30);

g6:=f6 div f2;
g6:=g6 div f3; 
g10:=f10 div f2;
g10:=g10 div f5; 

g15:=f15 div f3;
g15:=g15 div f5;

g30:=f30 div f2;
g30:=g30 div f3;
g30:=g30 div f5;
g30:=g30 div g6;
g30:=g30 div g10;
g30:=g30 div g15;

Factorisation(g30: DegreeLimit:=6);
//Output: []
//In this case g_30 does not have any irreducible factors of degree less then or equal to 6.


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

//E has a rational 9 isogeny and mod 2 Galois representation is surjective

E:=EllipticCurve([0,-27]);
E;
Q<x>:=PolynomialRing(Rationals());
K<w>:=NumberField(x^2+3);
E0:=BaseChange(E,K);
DescentInformation(E0);
//Torsion subgroup of E0 contains 12 points.
Points(E0:Bound:=30);
//This outputs 12 points on E0, so we've found them all.

//For each point (x,y) on E(K<w>), we have that if x is rational, then x is one 
// of the values 0,3 or -6, which is impossible.


//E does not have a rational 9-isogeny, has a rational 3-isogeny and G_E(2)=2B

A<x,y>:=AffineSpace(Rationals(),2);
C:=Curve(A,(x+3)*(x^2-3*x+9)*(x^3+3)^3*y-256*(y+1)^3*x^3);
D:=ProjectiveClosure(C);
R<x,y>:=PolynomialRing(Rationals(),2);
tr,x:=IsHyperelliptic(D);
x1,f:=SimplifiedModel(x);
f2:=Inverse(f);
x1;
J:=Jacobian(x1);
RankBound(J);

tr,g:=IsIsomorphic(D, x);
tr,g2:=IsInvertible(g);
g2:=Inverse(g);

pts:=Chabauty0(J);
for i:=1 to #pts do
g2(f2(pts[i]));
end for;

//The only points on D for which xy=/=0 are (3,-16,1) and (-3,-1,1). Plugging in y=-16 and y=-1 in j(E)=256(y+1)^3/y we get that 
// j(E)=54000 or j(E)=0. In both of these cases, E has CM.


_<t>:=RationalFunctionField(Rationals());
f1:=2^(10)*3^(3)*t^3*(1-4*t^3);  //Q(E[2])<=Q(E[3])
f2:=(t+3)*(t^2-3*t+9)*(t^3+3)^3/t^3; //E does not have a rational 9-isogeny, has a rational 3-isogeny
R<x,y>:=PolynomialRing(Rationals(),2);
C:=ProjectiveClosure(Curve(AffineSpace(R),Numerator(Evaluate(f1,x)-Evaluate(f2,y))));
Genus(C);
bound:=10^3;
Pts := PointSearch(C,bound);
Pts:=[p : p in Pts | Multiplicity(p) eq 1];
assert #Pts ne 0;
Pt:=Pts[1];
E,mp1 := EllipticCurve(C,Pt);
CremonaReference(E);
EE,mm := MinimalModel(E);
mm:=mm^-1;
MW:=AbelianInvariants(MordellWeilGroup(EE));
//printf "Abelian Invariants of  MW %o\n", MW;
DescentInformation(E);


T,mp2 := TorsionSubgroup(EE);
PtsC := { };
for p in T do
	PtsC := PtsC join RationalPoints(mm(mp2(p)) @@ mp1);
end for;
PtsC;
j:=f1;
{Evaluate(j,P[1]/P[3]) : P in PtsC | Evaluate(Denominator(j),P[1]) ne 0 and P[3] ne 0};

//Output: {0}. But since j(E)=0 implies that E has CM, we are done.

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

//C6xC12

//X_(3B)xX_20
_<t>:=RationalFunctionField(Rationals());
f1:=27*(t+1)*(t+9)^3/t^3; //3B.1.1 or 3B.1.2
f2:=(t^2-3)^3*(-4*t^2+32*t+44)/(t+1)^4; //X_20
R<x,y>:=PolynomialRing(Rationals(),2);
C:=ProjectiveClosure(Curve(AffineSpace(R),Numerator(Evaluate(f1,x)-Evaluate(f2,y))));

bound:=10^3;
Pts := PointSearch(C,bound);
Pts:=[p : p in Pts | Multiplicity(p) eq 1];
assert #Pts ne 0;
Pt:=Pts[1];
E,mp1 := EllipticCurve(C,Pt);
CremonaReference(E);
EE,mm := MinimalModel(E);
mm:=mm^-1;
MW:=AbelianInvariants(MordellWeilGroup(EE));
//printf "Abelian Invariants of  MW %o\n", MW;
DescentInformation(E);


T,mp2 := TorsionSubgroup(EE);
PtsC := { };
for p in T do
	PtsC := PtsC join RationalPoints(mm(mp2(p)) @@ mp1);
end for;
PtsC;
j:=f1;
{Evaluate(j,P[1]/P[3]) : P in PtsC | Evaluate(Denominator(j),P[1]) ne 0 and P[3] ne 0};
 	
//Output: {-35937/4, 109503/64}


//But Q(E[2]) is contained in Q(E[3])), so j(E)=2^10*3^3*t^3*(1-4*t^3) for some rational number t and 
// none of the two invariants listed above is of that form, which we now check.

P<t>:=PolynomialRing(Rationals());
f:=2^10*3^3*t^3*(1-4*t^3)-35937/4;
g:=2^10*3^3*t^3*(1-4*t^3)-109503/64;
Roots(f), Roots(g);
//Output gives us that no rational roots of these equation exist.


/////// X_(3Ns) x X_(20)
A<x,t>:=AffineSpace(Rationals(),2);
C:=Curve(A,x^3*(t+1)^4-(-4*t^2+32*t+44)*(t^2-3)^3);
D:=ProjectiveClosure(C);
R<x,y>:=PolynomialRing(Rationals(),2);
tr,x:=IsHyperelliptic(D);
x1,f:=SimplifiedModel(x);
f2:=Inverse(f);
x1;
J:=Jacobian(x1);
RankBound(J);

tr,g:=IsIsomorphic(D, x);
tr,g2:=IsInvertible(g);
g2:=Inverse(g);

pts:=Chabauty0(J);
for i:=1 to #pts do
g2(f2(pts[i]));
end for;
//////


//3Cs1.1 i 2Cn
_<t>:=RationalFunctionField(Rationals());
f1:=t^3; //3Cs.1.1
f2:=t^2+1728; //2Cn
R<x,y>:=PolynomialRing(Rationals(),2);
C:=ProjectiveClosure(Curve(AffineSpace(R),Numerator(Evaluate(f1,x)-Evaluate(f2,y))));

bound:=10^3;
Pts := PointSearch(C,bound);
Pts;
Pts:=[p : p in Pts | Multiplicity(p) eq 1];
assert #Pts ne 0;
Pt:=Pts[1];
E,mp1 := EllipticCurve(C,Pt);
CremonaReference(E);
EE,mm := MinimalModel(E);
mm:=mm^-1;
MW:=AbelianInvariants(MordellWeilGroup(EE));
//printf "Abelian Invariants of  MW %o\n", MW;
DescentInformation(E);


T,mp2 := TorsionSubgroup(EE);
PtsC := { };
for p in T do
	PtsC := PtsC join RationalPoints(mm(mp2(p)) @@ mp1);
end for;
j:=f1;
{Evaluate(j,P[1]/P[3]) : P in PtsC | Evaluate(Denominator(j),P[1]) ne 0 and P[3] ne 0};
//Output: {1728}


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

Z9:=Integers(9);
Z3:=Integers(3);
Sub9:=[H`subgroup: H in Subgroups(GL(2,Z9))];
Borel9:=sub<GL(2,Z9)|[1,1,0,1],[1,0,0,2],[8,0,0,8],[2,0,0,1]>; 


B311:=sub<GL(2,Z3)|[1,0,0,2],[1,0,1,1]>;
B312:=sub<GL(2,Z3)|[2,0,0,1],[1,0,1,1]>; 
Cs311:=sub<GL(2,Z3)|[1,0,0,2]>;

ImmB311:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,B311)];
ImmB311:=[H: H in ImmB311 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
ImmB311:=[H: H in ImmB311 | Order(H) le 6];
ImmB311;
//Output: []

ImmB312:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,B312)];
ImmB312:=[H: H in ImmB312 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
ImmB312:=[H: H in ImmB312 | Order(H) le 6];
ImmB312;
//Output: []

ImmCs311:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,Cs311)];
ImmCs311:=[H: H in ImmCs311 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
ImmCs311:=[H: H in ImmCs311| Order(H) le 6];
ImmCs311;
//Output: Three groups are given in the output, but they're all conjugate subgroups of Borel subgroup of GL(2,9). We check this.

for h in ImmCs311 do
	if IsConjugateSubgroup(GL(2,Z9),Borel9,h) eq false then
		h;
	end if;
end for;

//Previous loop outputs nothing, as expected.


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

//C25

Z25:=Integers(25);
Z5:=Integers(5);
Sub25:=[H`subgroup: H in Subgroups(GL(2,Z25))];

// Let G be a subgroup of GL(2,Z/nZ) acts on the left: M*v (M in G) 
// then we need to transpose to work in Magma 
// (since Magma the subgroups of GL(2,Z/nZ) acts on the right: v*M (M in G) 

G5B11:=sub<GL(2,Z5) | {[1,0,0,2],[1,0,1,1]}>; // Transpose the generators on Sutherland
G5Cs11:=sub<GL(2,Z5) | {[1,0,0,2]}>;
G5Cs13:=sub<GL(2,Z5)|{[3,0,0,4]}>;
G5Cs41:=sub<GL(2,Z5)|{[4,0,0,4],[1,0,0,2]}>;
G5B14:=sub<GL(2,Z5)|[[4,0,0,3],[1,0,1,1]]>;
G5B41:=sub<GL(2,Z5)|[[4,0,0,4],[1,0,1,1],[1,0,0,2]]>;


// Imm5B11 is the set of subgroups GG of GL(2,Z/25Z) (up to conjugacy) such that GG = G (mod 5), where G=5B.1.1
Im_rho:=G5B11;
Imm5B11:=[H  : H in Sub25 | IsConjugate(GL(2,Z5),sub<GL(2,Z5) | {GL(2,Z5)!m : m in Generators(H)}>,Im_rho)];

// Imm5Cs11 is the set of subgroups GG of GL(2,Z/25Z) (up to conjugacy) such that GG = G (mod 5), where G=5Cs.1.1
Im_rho:=G5Cs11;
Imm5Cs11:=[H  : H in Sub25 | IsConjugate(GL(2,Z5),sub<GL(2,Z5) | {GL(2,Z5)!m : m in Generators(H)}>,Im_rho)];

// Imm5Cs13 is the set of subgroups GG of GL(2,Z/25Z) (up to conjugacy) such that GG = G (mod 5), where G=5Cs.1.3
Im_rho:=G5Cs13;
Imm5Cs13:=[H  : H in Sub25 | IsConjugate(GL(2,Z5),sub<GL(2,Z5) | {GL(2,Z5)!m : m in Generators(H)}>,Im_rho)];

// Imm5Cs41 is the set of subgroups GG of GL(2,Z/25Z) (up to conjugacy) such that GG = G (mod 5), where G=5Cs.4.1
Im_rho:=G5Cs41;
Imm5Cs41:=[H  : H in Sub25 | IsConjugate(GL(2,Z5),sub<GL(2,Z5) | {GL(2,Z5)!m : m in Generators(H)}>,Im_rho)];

// Imm5B14 is the set of subgroups GG of GL(2,Z/25Z) (up to conjugacy) such that GG = G (mod 5), where G=5B.1.4
Im_rho:=G5B14;
Imm5B14:=[H  : H in Sub25 | IsConjugate(GL(2,Z5),sub<GL(2,Z5) | {GL(2,Z5)!m : m in Generators(H)}>,Im_rho)];

// Imm5B41 is the set of subgroups GG of GL(2,Z/25Z) (up to conjugacy) such that GG = G (mod 5), where G=5B.4.1
Im_rho:=G5B41;
Imm5B41:=[H  : H in Sub25 | IsConjugate(GL(2,Z5),sub<GL(2,Z5) | {GL(2,Z5)!m : m in Generators(H)}>,Im_rho)];




Imm:=Imm5B11;

G_25:={};
for GG in Imm do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],5),IsDivisibleBy(Eltseq(v)[2],5)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_25:=G_25 join names; 
end for;
G_25;


Imm:=Imm5Cs11;

G_25:={};
for GG in Imm do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],5),IsDivisibleBy(Eltseq(v)[2],5)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_25:=G_25 join names; 
end for;
G_25;

Imm:=Imm5B41;

G_25:={};
for GG in Imm do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],5),IsDivisibleBy(Eltseq(v)[2],5)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_25:=G_25 join names; 
end for;
G_25;

Imm:=Imm5B14;

G_25:={};
for GG in Imm do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],5),IsDivisibleBy(Eltseq(v)[2],5)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_25:=G_25 join names; 
end for;
G_25;

Imm:=Imm5Cs41;

G_25:={};
for GG in Imm do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],5),IsDivisibleBy(Eltseq(v)[2],5)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_25:=G_25 join names; 
end for;
G_25;

Imm:=Imm5Cs13;

G_25:={};
for GG in Imm do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],5),IsDivisibleBy(Eltseq(v)[2],5)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_25:=G_25 join names; 
end for;
G_25;


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

//C36

Z9:=Integers(9);
Z3:=Integers(3);
Sub9:=[H`subgroup: H in Subgroups(GL(2,Z9))];

B311:=sub<GL(2,Z3)|[1,0,0,2],[1,0,1,1]>; 
B312:=sub<GL(2,Z3)|[2,0,0,1],[1,0,1,1]>; 
Cs311:=sub<GL(2,Z3)|[1,0,0,2]>;
Cs3:=sub<GL(2,Z3)|[1,0,0,2],[2,0,0,2]>;
B3:=sub<GL(2,Z3)|[1,0,0,2],[2,0,0,2],[1,0,1,1]>;

ImmB311:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,B311)];
ImmB311:=[H: H in ImmB311 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
#ImmB311; //Ni jedna grupa nije reda 6

ImmB312:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,B312)];
ImmB312:=[H: H in ImmB312 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
#ImmB312; //Ni jedna grupa nije reda 6

ImmCs311:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,Cs311)];
ImmCs311:=[H: H in ImmCs311 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
#ImmCs311; //Ni jedna grupa nije reda 6

ImmCs3:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,Cs3)];
ImmCs3:=[H: H in ImmCs3 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
#ImmCs3;

ImmB3:=[H: H in Sub9 | IsConjugate(GL(2,Z3),sub<GL(2,Z3) | {GL(2,Z3)!m : m in Generators(H)}>,B3)];
ImmB3:=[H: H in ImmB3 | #{Determinant(g): g in sub<GL(2,Z9)|H,[8,0,0,8]>} eq 6];
#ImmB3;

G_9:={};
for GG in ImmB311 do
    V25:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],3),IsDivisibleBy(Eltseq(v)[2],3)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_9:=G_9 join names; 
end for;
G_9;

G_9:={};
for GG in ImmB312 do
    V9:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],3),IsDivisibleBy(Eltseq(v)[2],3)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V25 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_9:=G_9 join names; 
end for;
G_9;


G_9:={};
for GG in ImmCs311 do
    V9:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],3),IsDivisibleBy(Eltseq(v)[2],3)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V9 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6}; 
    G_9:=G_9 join names;
end for;
G_9;


G_9:={};
for GG in ImmCs3 do
    V9:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],3),IsDivisibleBy(Eltseq(v)[2],3)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V9 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_9:=G_9 join names; 
end for;
G_9;


G_9:={};
for GG in ImmB3 do
    V9:={ v : v in RSpace(GG) | not IsZero(v) and not &and[IsDivisibleBy(Eltseq(v)[1],3),IsDivisibleBy(Eltseq(v)[2],3)] };
    names:={GroupName(quo<GG | Core(GG,Stabiliser(GG,v))>) : v in V9 | Integers()!(Order(GG)/Order(Stabiliser(GG,v))) eq 6};
    G_9:=G_9 join names; 
end for;
G_9;

//All the groups are of generalized S3-type. It follows that Q(P_9)<=Q(3^inf).


for rzb in RZB do
	T:=rzb[3];
	if T[1][2] eq 6 then
	rzb[1];
	end if;
end for;

//Output: X20,X20a. Hence, Q(P4) is S3 extension of Q, so Q(P_4)<=Q(3^inf).

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

//C45

// If E/Q has a rational 15-isogeny, then j(E)=a_i, for i<=4. Everything else is the same as in C63 case.
a1:=-5^2/2;
a2:=-5^2*241^3/2^3;
a3:=-29^3*5/2^(5);
a4:=211^3*5/2^(15);

list:=<a1,a2,a3,a4>;

for a in list do
E:=EllipticCurveWithjInvariant(a);

g3:=DivisionPolynomial(E,3);
g5:=DivisionPolynomial(E,5);
f9:=DivisionPolynomial(E,9);
f15:=DivisionPolynomial(E,15);

g9:=f9 div g3;

g15:=f15 div g5;
g15:=g15 div g3;

f45:=DivisionPolynomial(E,45);

g45:=f45 div g15;
g45:=f45 div g9;
g45:=f45 div g5;
g45:=f45 div g3;

//Now g45 is a polynomial whose roots are x-coordinates of point of exact order 45 on E.


Factorisation(g45: DegreeLimit:=6);

end for;

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

//C63

// If E/Q has a rational 21-isogeny, then j(E)=a_i, for i<=4.
a1:=-3^2*5^6/2^3;
a2:=3^3*5^3/2^1;
a3:=3^3*5^3*101^3/2^(21);
a4:=-3^3*5^3*383^3/2^7;

list:=<a1,a2,a3,a4>;

//take a random elliptic curve with j(E)=a
for a in list do
E:=EllipticCurveWithjInvariant(a); 

//computing  primitive division polynomial g63.
g3:=DivisionPolynomial(E,3);
g7:=DivisionPolynomial(E,7);
f9:=DivisionPolynomial(E,9);
f21:=DivisionPolynomial(E,21);

g9:=f9 div g3; //this is the 9th primitive division polynomial

g21:=f21 div g7;
g21:=g21 div g3; //this is the 21st primitive division polynomial

f63:=DivisionPolynomial(E,63);

g63:=f63 div g21;
g63:=g63 div g9;
g63:=g63 div g7;
g63:=g63 div g3; //this is the 63rd primitive division polynomial

//g63 is a polynomial whose roots are x-coordinates of point of exact order 63 on E.


Factorisation(g63: DegreeLimit:=6);

end for;



//C42
//We do exactly the same as in the C63 case.

a1:=-3^2*5^6/2^3;
a2:=3^3*5^3/2^1;
a3:=3^3*5^3*101^3/2^(21);
a4:=-3^3*5^3*383^3/2^7;

list:=<a1,a2,a3,a4>;

for a in list do
E:=EllipticCurveWithjInvariant(a);

g2:=DivisionPolynomial(E,2);
g3:=DivisionPolynomial(E,3);
g7:=DivisionPolynomial(E,7);
f6:=DivisionPolynomial(E,6);
f21:=DivisionPolynomial(E,21);

g6:=f6 div g2;
g6:=g6 div g3; 

g14:=f14 div g2;
g14:=g14 div g7;

g21:=f21 div g7;
g21:=g21 div g3;

f42:=DivisionPolynomial(E,42);

g42:=f42 div g21;
g42:=g42 div g14;
g42:=g42 div g7;
g42:=g42 div g6;
g42:=g42 div g3;
g42:=g42 div g2;

//Now g42 is a polynomial whose roots are x-coordinates of point of exact order 42 on E.


Factorisation(g42: DegreeLimit:=6);

end for;