/* p=3, 19-torsion
The following code proves that there is no 19-torsion over Q_\infty_3 */

E := EllipticCurve("361a2");
fac:=Factorization(DivisionPolynomial(E,19));
K:=CyclotomicField(27);
L:=Subfields(K,3)[1][1];
E1:=ChangeRing(E,L);
fac:=Factorization(DivisionPolynomial(E1,19));
fac;


/* p=3, 14-torsion 
 The following code proves that there is no 14-torsion over Q_\infty_3 */

E := EllipticCurve([1,0,1,-1,0]);
K:=CyclotomicField(9);
L:=Subfields(K,3)[1][1];
E1:=ChangeRing(E,L);
TorsionSubgroup(E1);
DescentInformation(E1);


/* p=3, 27-torsion 
 The following code finds the unique curve with 27-torsion over Q_\infty,3 */

E := EllipticCurve("27a2");
K:=CyclotomicField(27);
L:=Subfields(K,9)[1][1];
E1:=ChangeRing(E,L);
fac:=Factorization(DivisionPolynomial(E1,27) div DivisionPolynomial(E1,9));
fac;
w:=Roots(fac[1][1])[1][1];
w;

/* w:=6*L.1^7 - 93*L.1^6 + 576*L.1^5 - 1809*L.1^4 + 3015*L.1^3 - 2547*L.1^2 + 954*L.1
    - 123 */
_<y>:=PolynomialRing(L);
p:=Factorization(Discriminant(y^2 + y -(w^3 - 270*w - 1708))*RingOfIntegers(L));
p1:=p[1][1];
Norm(p1);
TorsionSubgroup(QuadraticTwist(E1,-3));



/*p=3, 18-torsion*/

/*This shows that the rank of J_1(18) over Q_1,3 is 0 */
_<x>:=PolynomialRing(Rationals());
C:=HyperellipticCurve(x^6+2*x^5+5*x^4+10*x^3+10*x^2+4*x+1);
K:=CyclotomicField(9);
L:=Subfields(K,3)[1][1];
C1:=ChangeRing(C,L);
J:=Jacobian(C1);
RankBound(J);

/*This shows that the torsion J_1(18) over Q_1,3 is contained in  Z/21Z x Z/21Z */

O:=RingOfIntegers(L);
d:=17;
Factorization(d*O);
Cr:=ChangeRing(C,GF(17));
Jr:=Jacobian(Cr);
AbelianGroup(Jr);

/*This shows that the torsion J_1(18) over Q_1,3 is contains in  Z/7Z x Z/21Z */

pts:=Points(C1:Bound:=100);
S:={};
A:=[];
P1:=pts[4]-pts[7];
P2:=pts[8]-pts[9];
for i:=1 to 21 do 
for j:=1 to 21 do 
S:=S join {i*P1+j*P2};
end for;
end for;
S;

/*This shows that the only point on X_1(18) over Q_1,3 are the 12 cusps */

for a in S do p:=Roots(a[1]); 
p;
if IsEmpty(p) eq false then 
for i:=1 to #p do
A:=A cat [p[i][1]];
end for;
end if;
emd for;

P2:=PointsAtInfinity(C1);
for i := 1 to #A do P2:=P2 join Points(C1, A[i]); end for; 
assert P2 eq pts;


/*p=2, 17-torsion*/

E := EllipticCurve("14450p2");
_<x>:=PolynomialRing(Rationals());
fac:=Factorization(DivisionPolynomial(E,17));
p1:=fac[1][1];
p2:=fac[2][1];
L:=CyclotomicField(32);
a:=Subfields(L,8);
K:=a[2][1];
IsTotallyReal(K);
IsCyclcic(GaloisGroup(K));
P<y>:=PolynomialRing(K);
pz:=ChangeRing(p1,P);
Factorization(pz);
pz:=ChangeRing(p2,P);
Factorization(pz);
Factorization(y^4 + 3372*y^3 + 1316094*y^2 - 836445572*y - 267892141399);
Factorization(y^4 + 482*y^3 + 1144*y^2 - 15809842*y - 958623689);
K:=a[3][1];
IsTotallyReal(K);
IsCyclic(GaloisGroup(K));
P<y>:=PolynomialRing(K);
Factorization(y^4 + 3372*y^3 + 1316094*y^2 - 836445572*y - 267892141399);
Factorization(y^4 + 482*y^3 + 1144*y^2 - 15809842*y - 958623689);


E := EllipticCurve("14450p1");
fac:=Factorization(DivisionPolynomial(E,17));
p1:=fac[1][1];
L:=CyclotomicField(32);
a:=Subfields(L,8);
K:=a[2][1];
IsTotallyReal(K);
IsCyclic(GaloisGroup(K));
P<x>:=PolynomialRing(K);
pz:=ChangeRing(p1,P);
Factorization(pz);
Factorization(x^8 - 226*x^7 + 18372*x^6 - 543828*x^5 - 9242705*x^4 + 1127218758*x^3 -
    33006143963*x^2 + 437271444481*x - 2252576338909);
K:=a[3][1];
IsTotallyReal(K);
IsCyclic(GaloisGroup(K));
P<x>:=PolynomialRing(K);
Factorization(x^8 - 226*x^7 + 18372*x^6 - 543828*x^5 - 9242705*x^4 + 1127218758*x^3 -
    33006143963*x^2 + 437271444481*x - 2252576338909);


/* p=2, 15-torsion */

E := EllipticCurve([1, 1, 1, 0, 0]);
DescentInformation(E);
_<x>:=PolynomialRing(Rationals());
K:=NumberField(x^4-4*x^2+2);
E1:=ChangeRing(E,K);
DescentInformation(E1);

/* p=2, 2x10-torsion */

E := EllipticCurve([0, 1, 0, -1, 0]);
DescentInformation(E);
_<x>:=PolynomialRing(Rationals());
K:=NumberField(x^2-2);
E1:=ChangeRing(E,K);
DescentInformation(E1);





/*Code for Theorem 6.1 */

/* The 1-parameter families for curves with a point of order 4,6,8 respectively
taken from Lozano-Robledo, "Elliptic Curves, Modular Forms, and Their L-functions." */
F<t> := FunctionField(Rationals());
a4 := 0;
b4 := t;
a6 := t;
b6 := t^2+t;
a8 := (2*t-1)*(t-1)/t;
b8 := (2*t-1)*(t-1);

/* Define the generic curves with a point of order n over Q. */
E4 := EllipticCurve([1-a4,-b4,-b4,0,0]); //has a rational point of order 4
E6 := EllipticCurve([1-a6,-b6,-b6,0,0]); //has a rational point of order 6
E8 := EllipticCurve([1-a8,-b8,-b8,0,0]); //has a rational point of order 8

/*Check the torsion subgroup is what we think it should be */
assert IsIsomorphic(TorsionSubgroup(E4),AbelianGroup([4]));
assert IsIsomorphic(TorsionSubgroup(E6),AbelianGroup([6]));
assert IsIsomorphic(TorsionSubgroup(E8),AbelianGroup([8]));

/* Compute the discriminant mod squares */
assert IsSquare((16*t+1)*Discriminant(E4)); 
assert IsSquare((9*t+1)*(t+1)*Discriminant(E6));
assert IsSquare(2*(t^2-t+1/8)*Discriminant(E8));

/* Define the curves whose t coordinates would give specialization with discriminant 2 mod 
squares. */
A<t,s> := AffineSpace(Rationals(),2);
C42 := ProjectiveClosure(Curve(A,(16*t+1) - 2*s^2));
C62 := ProjectiveClosure(Curve(A,(9*t+1)*(t+1)- 2*s^2));
C82 := ProjectiveClosure(Curve(A,2*(t^2-t+1/8) - 2*s^2));

/* Define the curves whose t coordinates would give specialization with discriminant -1 mod 
squares. */
C41 := ProjectiveClosure(Curve(A,(16*t+1) + s^2));
C61 := ProjectiveClosure(Curve(A,(9*t+1)*(t+1) + s^2));
C81 := ProjectiveClosure(Curve(A,2*(t^2-t+1/8) + s^2));

//Show these are genus zero curves.
assert Genus(C42) eq 0;
assert Genus(C62) eq 0;
assert Genus(C82) eq 0;
assert Genus(C41) eq 0;
assert Genus(C61) eq 0;
assert Genus(C81) eq 0;

//Check that each of these curves has points on it.
assert #Points(C42 : Bound := 10^3) ge 1;
assert #Points(C62 : Bound := 10^3) ge 1;
assert #Points(C82 : Bound := 10^3) ge 1;
assert #Points(C41 : Bound := 10^3) gt 1;
assert #Points(C61 : Bound := 10^3) gt 1;
assert #Points(C81 : Bound := 10^3) gt 1;

/* Defined the family of elliptic curves with a point of order 2 and check that the torsion 
group is what we say it is. */
F<t> := FunctionField(Rationals());
E2 := EllipticCurve([0,4*(1-2*t^2)^-1,0,4*(1-2*t^2)^-1,0]); 
assert IsIsomorphic(TorsionSubgroup(E2),AbelianGroup([2]));

//Show the discriminant is 2 mod squares. 
assert IsSquare(2*Discriminant(E2));

//Show the family is non-isotrivial.
assert Degree(Numerator(jInvariant(E2))) ge 0;




















