// We first need the groups and rational functions of the paper to be loaded.
load "Equations.txt";

function ApplicableSubgroups(ell, G)
        // Input: a prime ell and a subgroup G of GL_2(F_ell).
	// Output: all proper applicable subgroups of G up to conjugacy in GL_2(F_ell).
    	M:=[H`subgroup: H in Subgroups(G)];
        S:={{H1: H1 in M | IsConjugate(GL(2,ell),H1,H2)}: H2 in M};
        M:=[Random(s): s in S];   // list of subgroups of G up to conjugacy in GL_2(F_ell)
    	M:=[H: H in M | #{Determinant(h): h in H} eq ell-1 and GL(2,ell)![-1,0,0,-1] in H and H ne G];
    	return [H: H in M |  #{h: h in H | Determinant(h) eq -1 and Trace(h) eq 0 } ne 0];
end function;

function AllApplicableSubgroups(ell,G)
	// Determines if a sequence G of subgroups of GL_2(F_ell) are representatives of
	// the distinct conjugacy classes of applicable subgroups of GL_2(F_ell).
	H:=ApplicableSubgroups(ell, GL(2,ell));
	if #H ne #G then return false; end if;
	return &and[&or{IsConjugate(GL(2,ell),G[i],H[j]): j in [1..#H]} : i in [1..#G]];
end function;

function FindRelation(h,j,n,K)
	// Input: Puiseux series h and j with coefficients in K such that [K(h):K(j)]=n.
	// Output: J(t) in K(t) such that j=J(h).
	// Will produce an error if not enough terms of j and h are known.
	P<[a]>:=PolynomialRing(K,2*(n+1));
	R<q>:=PuiseuxSeriesRing(P);
        h:=R!h; j:=R!j;

	s:= j*&+[a[i+1]*h^i : i in [0..n]] - &+[a[i+n+2]*h^i : i in [0..n]];
	A:=[];
	repeat
	    v:=Valuation(s); c:=P!Coefficient(s,v);   s:=s-c*q^v;
	    A:= A cat [ [K!MonomialCoefficient(c,a[i]) : i in [1..2*n+2]] ];
	    B:=Matrix(K,A);  //Rank(B);
        until Rank(B) eq 2*n+1;
        v:=Basis(NullSpace(Transpose(B)))[1];

        F<t>:=FunctionField(K);
        return (&+[v[i+n+2]*t^i : i in [0..n]])/(&+[v[i+1]*t^i : i in [0..n]]);
end function;


// ell=2
print "Performing ell=2 verifications.";
assert AllApplicableSubgroups(2,G2)  and  [Index(GL(2,2),H): H in G2] eq [Degree(j): j in J2];
_<q>:=PuiseuxSeriesRing(Rationals():Precision:=100);
h1:=16*DedekindEta(q^2)^8/DedekindEta(q^(1/2))^8;
assert FindRelation(h1,jInvariant(q),6,Rationals()) eq J2[1];
assert Evaluate(J2[2],t^2/(t+1)) eq J2[1];
assert Evaluate(J2[3],(16*t^3 + 24*t^2 - 24*t - 16)/(t^2+t)) eq J2[1];

// ell=3
print "Performing ell=3 verifications.";
assert AllApplicableSubgroups(3,G3);
assert [Index(GL(2,3),H): H in G3] eq [12,6,4,3] and [Degree(j): j in J3] eq [12,6,4,3];
_<q>:=PuiseuxSeriesRing(Rationals():Precision:=100);
h1:= 1/3*DedekindEta(q^(1/3))^3/DedekindEta(q^3)^3;
assert FindRelation(h1,jInvariant(q),12,Rationals()) eq J3[1];
assert Evaluate(J3[2],(t^2+3*t+3)/t) eq J3[1];
assert Evaluate(J3[3],t*(t^2+3*t+3)) eq J3[1];
assert Evaluate(J3[4],3*(t+1)*(t-3)/t) eq J3[2];

// ell=5
print "Performing ell=5 verifications.";
assert AllApplicableSubgroups(5,G5);
assert [Index(GL(2,5),H): H in G5] eq [60,30,30,15,12,12,10,6,5] and [Degree(j): j in J5] eq [60,30,30,15,12,12,10,6,5];

r:=q^(1/5)*&*[ (1-q^i)^KroneckerSymbol(5,i) : i in [1..30] ]+O(q^30);  h1:=1/r;
assert FindRelation(h1,jInvariant(q),60,Rationals()) eq J5[1];

K<d>:=QuadraticField(5); w:=(1+d)/2;
F<t>:=FunctionField(K);
assert Evaluate(J5[2],t-1-1/t) eq J5[1];
assert Evaluate(J5[3],((-3+w)*t-5)/(t+3-w)) eq J5[2];
assert Evaluate(J5[4],t+5/t) eq J5[2];
assert Evaluate(J5[5],t^5) eq J5[1];
assert Evaluate(J5[6],(-(w-1)^5*t+1)/(t+(w-1)^5)) eq J5[5];
assert Evaluate(J5[7],-(t^3+10*t^2+25*t+25)/(2*t^3+10*t^2+25*t+25)) eq J5[3];
assert Evaluate(J5[8],(t-11-1/t)/25) eq J5[5];
assert Evaluate(J5[9],(t+5)*(t^2-5)/(t^2+5*t+5)) eq J5[4];

// ell=7
print "Performing ell=7 verifications.";
G7extra:=[
	sub<GL(2,7) | {[-1,0,0,-1],[1,0,0,3]}>,
	sub<GL(2,7) | {[1,0,0,-1],[3,0,0,3]}>,
	sub<GL(2,7) | {[0,-2,2,0],[1,0,0,-1]}>,        
	sub<GL(2,7) | {[1,0,0,3],[3,0,0,1]}>,
        sub<GL(2,7) | {[1,-1,1,1],[1,0,0,-1]}>
	];
assert AllApplicableSubgroups(7,G7 cat G7extra);
assert [Index(GL(2,7),H): H in G7 cat G7extra] eq [56,28,24,24,24,21,8,168,168,84,56,42];

F<t>:=FunctionField(Rationals());
M:=50;  R<q>:=PuiseuxSeriesRing(Rationals():Precision:=M);
x:=-q^(4/7)*&*[(1-q^n)^3*(1-q^(7*n)) +O(q^M): n in [1..M]]*&*[1-q^n +O(q^M): n in [1..M] | n mod 7 in {1,6}];
y:= q^(2/7)*&*[(1-q^n)^3*(1-q^(7*n)) +O(q^M): n in [1..M]]*&*[1-q^n +O(q^M): n in [1..M] | n mod 7 in {2,5}];
z:= q^(1/7)*&*[(1-q^n)^3*(1-q^(7*n)) +O(q^M): n in [1..M]]*&*[1-q^n +O(q^M): n in [1..M] | n mod 7 in {4,3}];

h4:=-y^2*z/x^3; 
F1:=t+1/(1-t)+(t-1)/t-8;
h7:=Evaluate(F1,h4);
assert FindRelation(h7,jInvariant(q),8,Rationals()) eq J7[7];
assert Evaluate(J7[7],F1) eq J7[4];
K<z7>:=CyclotomicField(7);
F<t>:=FunctionField(K);
beta:=4+3*z7+3/z7+z7^2+1/z7^2;  gamma:=z7^5+z7^4+z7^3+z7^2+1;
F2:=(beta*t-(beta-1))/(t-beta);   F3:=(t-gamma)/(gamma*t-(gamma-1));
assert Evaluate(J7[3],F2) eq J7[4];
assert Evaluate(J7[5],F3) eq J7[4];

_<x,y,z>:= PolynomialRing(Rationals(),3);
H7:=Evaluate(F1,-y^2*z/x^3);
H2:=(x+y+z)^2/(x*y+y*z+z*x);
V:=(x*y+y*z+z*x)^2/((x+y+z)*x*y*z);
f:=Numerator(H7+(H2^3-4*H2^2+3*H2+1)*((H2^2-5*H2+1)*V+7));
assert Denominator(f/(x^3*y+y^3*z+z^3*x)) eq 1;

FF<h>:=FunctionField(Rationals());  RR<w>:=PolynomialRing(FF);
_<w>:=quo<RR| w^2-(h^4-10*h^3+27*h^2-10*h-27)>;
hh:=1/2*(h^3-4*h^2+3*h+1)*((h^4-10*h^3+27*h^2-10*h-13)-(h^2-5*h+1)*w);
assert Evaluate(J7[7],hh) eq Evaluate(J7[2],h);

P2<x,y,z>:=ProjectiveSpace(Rationals(),2);
E1:=EllipticCurve([0,0,0,-1715,33614]); E2:=EllipticCurve([0,0,0,-29155,1915998]);
assert #MordellWeilGroup(E1) eq 2 and #MordellWeilGroup(E2) eq 2;
C:=Curve(P2,y^2*z^2-(x^4-10*x^3*z+27*x^2*z^2-10*x*z^3-27*z^4));
assert IsIsomorphic(EllipticCurve(C,C![0,1,0]), QuadraticTwist(E1,-7));
C:=Curve(P2,-7*y^2*z^2-(x^4-10*x^3*z+27*x^2*z^2-10*x*z^3-27*z^4));
assert IsIsomorphic(E1,EllipticCurve(C,C![5/2,1/4,1]));
C:=Curve(P2,y^2*z^2-7*(2*x^4-14*x^3*z+21*x^2*z^2+28*x*z^3+7*z^4));
assert IsIsomorphic(E2,EllipticCurve(C,C![0,7,1]));

// ell=11
print "Performing ell=11 verifications.";
Pol<x>:=PolynomialRing(Rationals());
E:=EllipticCurve([1,1,1,-305,7888]);
assert jInvariant(E) eq -11^2 and Conductor(E) eq 11^2;
assert TraceOfFrobenius(E,2) eq -1;
f:=x^5-129*x^4+800*x^3+81847*x^2-421871*x-4132831;
F:=Factorization(DivisionPolynomial(E,11));
assert #F eq 2 and Degree(F[2][1]) eq 55 and F[1][1] eq f;
K<w>:=quo<Pol|f>;
E:=BaseChange(E,K);
P:=E![w,-(w^4-79*w^3-3150*w^2+12193*w+1520110)/11^4];
assert 11*P eq E![0,1,0];

E:=EllipticCurve([1,1,0,-3632,82757]);
assert jInvariant(E) eq -11*131^3 and Conductor(E) eq 11^2;
assert TraceOfFrobenius(E,2) eq 1;
f:=x^5-129*x^4+4793*x^3+9973*x^2-3694800*x+52660939;
F:=Factorization(DivisionPolynomial(E,11));
assert #F eq 2 and Degree(F[2][1]) eq 55 and F[1][1] eq f;
K<w>:=quo<Pol|f>;
E:=BaseChange(E,K);
P:=E![w,-(w^4-79*w^3+843*w^2+45468*w-722625)/11^3];
assert 11*P eq E![0,1,0];

K<x>:=FunctionField(Rationals());
P<Y>:=PolynomialRing(K);
_<y>:=quo<P|Y^2+Y-(x^3-x^2-7*x+10)>;
f1:=x^2+3*x-6;
f2:=11*(x^2-5)*y+(2*x^4+23*x^3-72*x^2-28*x+127);
f3:=6*y+11*x-19;
f4:=22*(x-2)*y+(5*x^3+17*x^2-112*x+120);
f5:=11*y+(2*x^2+17*x-34);
f6:=(x-4)*y-(5*x-9);
J:=(f1*f2*f3*f4)^3/(f5^2*f6^11);
f:=P!((K!A11)*J);
a:=Coefficient(f,1);
b:=Coefficient(f,0); 
A:=(K!A11); B:=-2*b+a; C:=(b^2-b*a-a^2*(x^3-x^2-7*x+10))/A;
assert A*J eq a*y+b;
assert B eq K!B11;
assert C eq K!C11;
assert A*J^2+B*J+C eq 0;
Delta:=K!(B^2-4*A*C);
assert IsSquare(Delta*(x^3-x^2-7*x+41/4));
assert #Roots(PolynomialRing(Rationals())!Delta) eq 0;
assert #Roots(PolynomialRing(Rationals())!a) eq 0;


// ell=13
print "Performing ell=13 verifications.";
K<z>:=CyclotomicField(13);
FF<t>:=FunctionField(K);

F2:=(t^3 - 4*t^2 + t + 1)/(t^2 - t);
F5:=(t^2 - 3*t - 1)/t;
assert Evaluate(J13[6],F2) eq J13[2];
assert Evaluate(J13[6],F5) eq J13[5];

alpha:= -z^11-z^10-z^3-z^2+1;
F1:=13*(t^2-t)/(t^3-4*t^2+t+1); phi1:=(alpha*t+(1-alpha))/(t-alpha);
assert Evaluate(F1,phi1) eq F2;
assert Evaluate(J13[6],F1) eq J13[1];

beta:=z^11+z^10+z^9+z^7+z^6+z^4+z^3+z^2+2;
F3:=(-5*t^3+7*t^2+8*t-5)/(t^3-4*t^2+t+1); phi3:= (beta*t-1)/(t+beta-1);
assert Evaluate(F3,phi3) eq F2;
assert Evaluate(J13[6],F3) eq J13[3];
gamma:=-z^11-z^8-z^7-z^6-z^5-z^2;
F4:=13*t/(t^2-3*t-1);  phi4:=((2-gamma)*t+1)/(t-2+gamma);
assert Evaluate(F4,phi4) eq F5;
assert Evaluate(J13[6],F4) eq J13[4];


//ell=17
print "Performing ell=17 verifications.";
E:=EllipticCurve([1,0,1,-190891,-36002922]);
assert jInvariant(E) eq -17*373^3/2^17 and Conductor(E) eq 2*5^2*17^2;
Pol<x>:=PolynomialRing(Rationals());
f:=x^4 + 482*x^3 + 1144*x^2 - 15809842*x - 958623689;
f2:=x^4 + 3372*x^3 + 1316094*x^2 - 836445572*x - 267892141399;
P:=DivisionPolynomial(E,17);
assert P mod (f*f2) eq 0;
assert IsIrreducible(P div (f*f2)) and Degree(P div (f*f2)) eq 8*17;
assert #GaloisGroup(f) eq 4;

Pol2<X>:=PolynomialRing(GF(17));
assert &join[{Order(r[1]): r in Roots(X^2-TraceOfFrobenius(E,p)*X+p)}: p in {103,137,307}] eq {4};

EE:=IsogenyFromKernel(E,f*f2);
assert IsIsomorphic(EE,EllipticCurve([1,0,1,-3041,64278]));
assert jInvariant(EE) eq -17^2*101^3/2;


//ell=37
print "Performing ell=37 verifications.";
E:=EllipticCurve([1,1,1,-8,6]);
assert jInvariant(E) eq -7*11^3 and Conductor(E) eq 5^2*7^2;
Pol<x>:=PolynomialRing(Rationals());
P:=DivisionPolynomial(E,37);
f:=x^6 - 15*x^5 - 90*x^4 - 50*x^3 + 225*x^2 + 125*x - 125;
f2:=x^6 - 15*x^5 + 50*x^4 + 125*x^3 - 825*x^2 + 1000*x - 125;
f3:=x^6 - 85*x^5 + 435*x^4 - 750*x^3 + 400*x^2 + 125*x - 125;
assert IsIrreducible(f) and IsIrreducible(f2) and IsIrreducible(f3);
assert P mod (f*f2*f3) eq 0;
P2:=P div (f*f2*f3);
assert IsIrreducible(P2) and Degree(P2) eq 18*37;
assert #GaloisGroup(f) eq 6;

EE:=IsogenyFromKernel(E,f*f2*f3);
assert IsIsomorphic(EE,EllipticCurve([1,1,1,-208083,-36621194]));
assert jInvariant(EE) eq -7*137^3*2083^3;

print "Done.";