/***************************************************************************************************
          PRINCIPALLY POLARIZED SQUARES OF ELLIPTIC CURVES WITH FIELD OF MODULI EQUAL TO Q
                    ALEXANDRE GÉLIN, EVERETT W. HOWE, AND CHRISTOPHE RITZENTHALER
***************************************************************************************************/

// This is the main part of our code, once the 1226 polarizations are computed.
// For each polarization, we check if there exist matrices P such that Q is field of moduli.
// The second step corresponds to the recovery of the rational invariants once it has been shown
// that Q is field of moduli. We also try to recover a model over Q when it is possible.

load "Pol.m"; //load the information about the polarizations

ListPol:=[**];
ListG2:=[];

//list of discriminants = 0 mod 4 classified by class number
List0:=[-4,-8];
List1:=[-20,-24,-40,-52,-88,-148,-232];
List2:=[-84,-120,-132,-168,-228,-280,-312,-340,-372,-408,-520,-532,-708,-760,-1012];
List3:=[-420,-660,-840,-1092,-1320,-1380,-1428,-1540,-1848];
List4:=[-5460];
ListAll:=List0 cat List1 cat List2 cat List3 cat List4;
//the ones that lead to a curve C
List:=[-8,-20,-24,-40,-52,-88,-148,-232,-84,-120,-132,-228,-312,-520,-708,-760];

//test if Q is field of moduli of E^2 or not
time for Disc in List do printf "\nDisc=%o. ",Disc;
    QF<u>:=QuadraticField(Disc); w:=u; FD:=-Disc/4;
    PolSet:=polar(Disc,w);
    printf "#Pol=%o.\n",#PolSet;
    G,f:=ClassGroup(QF); e:=Valuation(#G,2); 
    time for pol in PolSet do
	a:=IntegerRing()!pol[1]; d:=IntegerRing()!pol[2]; b:=pol[3];
	for k:=e to 1 by -1 do
	    I:=f(G.k); n:=Norm(I);
	    B1:=Floor(a*Sqrt(n)); B2:=Floor(a*Sqrt(n/FD)); N:=n*a^2;
	    NormSolutions:=[]; //solutions of the norm equation
	    for X1:=-B1 to B1 do
		for X2:=-B2 to B2 do
 		    for Z1:=-B1 to B1 do
			for Z2:=-B2 to B2 do
	 		    if X1^2+FD*X2^2+Z1^2+FD*Z2^2 eq N then
				Append(~NormSolutions,[(X1+w*X2)/a,(Z1+w*Z2)/a]);
			    end if;
			end for;
		    end for;
		end for;
	    end for;
	    ns:=#NormSolutions;
	    for i in [1..ns] do
		x:=NormSolutions[i][1]; z:=NormSolutions[i][2];
		for j in [i..ns] do
		    y:=NormSolutions[j][1]; t:=NormSolutions[j][2];
		    if Conjugate(x)*y+Conjugate(z)*t eq 0 and a*z in I and (x-b*z) in I
		       and (b*z+t) in I and 1/a*(b*x+y-b^2*z-b*t) in I then continue k; end if;
		end for;
	    end for;
	    printf("N "); continue pol; //Q is not field of moduli -> go to next polarization
	end for;
	printf("Y "); Append(~ListPol,<w,pol>); //Q is field of moduli -> store the pol
    end for;
end for;



//list of discriminants = 1 mod 4 classified by class number
List0:=[-3,-7,-11,-19,-43,-67,-163];
List1:=[-15,-35,-51,-91,-115,-123,-187,-235,-267,-403,-427];
List2:=[-195,-435,-483,-555,-595,-627,-715,-795,-1435];
List3:=[-1155,-1995,-3003,-3315];
ListAll:=List0 cat List1 cat List2 cat List3;
//the ones that lead to a curve C
List:=[-11,-19,-43,-67,-163,-595,-795];

//test if Q is field of moduli of E^2 or not
for Disc in List do printf "\nDisc=%o. ",Disc;
    QF<u>:=QuadraticField(Disc); w:=(1+u)/2; FD:=-Disc;
    PolSet:=polar(Disc,w);
    printf "#Pol=%o.\n",#PolSet;
    G,f:=ClassGroup(QF); e:=Valuation(#G,2); 
    time for pol in PolSet do
	a:=IntegerRing()!pol[1]; d:=IntegerRing()!pol[2]; b:=pol[3];
	for k:=e to 1 by -1 do
	    I:=f(G.k); n:=Norm(I);
	    B2:=Floor(2*a*Sqrt(n/FD)); N:=4*n*a^2;
	    NormSolutions:=[]; //solutions of the norm equation
	    for X2:=-B2 to B2 do
		for X1:=-Floor(a*Sqrt(n)-X2/2) to Floor(a*Sqrt(n)-X2/2) do
		    for Z2:=-B2 to B2 do
			for Z1:=-Floor(a*Sqrt(n)-Z2/2) to Floor(a*Sqrt(n)-Z2/2) do
			    if (2*X1+X2)^2+FD*X2^2+(2*Z1+Z2)^2+FD*Z2^2 eq N then
				Append(~NormSolutions,[(X1+w*X2)/a,(Z1+w*Z2)/a]);
			    end if;
			end for;
		    end for;
		end for;
	    end for;
	    ns:=#NormSolutions;
	    for i in [1..ns] do
		x:=NormSolutions[i][1]; z:=NormSolutions[i][2];
		for j in [i..ns] do
		    y:=NormSolutions[j][1]; t:=NormSolutions[j][2];
		    if Conjugate(x)*y+Conjugate(z)*t eq 0 and a*z in I and (x-b*z) in I
		       and (b*z+t) in I and 1/a*(b*x+y-b^2*z-b*t) in I then continue k; end if;
		end for;
	    end for;
	    printf("N "); continue pol; //Q is not field of moduli -> go to next polarization
	end for;
	printf("Y "); Append(~ListPol,<w,pol>); //Q is field of moduli -> store the pol
    end for;
end for;

//list all the 46 polarizations that correspond to a genus-2 curve
#ListPol;





//RECOVERY OF THE INVARIANTS
load "RedPerMat.m"; //load the tools for the period matrices
Prec:=220; //precision we need
SetDefaultRealField(RealField(Prec));
CF<I>:=ComplexField(Prec);
PR<X>:=PolynomialRing(RationalField());

time for Pol in ListPol do 
    w:=Pol[1]; pol:=Pol[2];
    M:=Matrix(Parent(w),2,2,[pol[1],pol[3],Conjugate(pol[3]),pol[2]]);
    Tau:=PerMat(M,w); //the period matrx
    TauRed:=Reduc(Tau); //the reduced one
    TauCF:=ChangeRing(TauRed,CF);
    S:=RosenhainInvariants(TauCF);
    f:=X*(X-1)*&*{X-a : a in S};
    C:=HyperellipticCurve(f);
    G2a:=[Real(g2) : g2 in G2Invariants(C)]; //approximated
    G2:=[BestApproximation(g,3^Prec) : g in G2a]; //recovered
    /* Factorisation(Denominator(G2[1])); */
    /* Factorisation(Denominator(G2[2])); */
    /* Factorisation(Denominator(G2[3])); */
    /* "\n"; */
    Append(~ListG2,G2);
end for;

//list all the 46 triplets of invariants of the genus-2 curves
#ListG2;

//for the equations of the 13 curves defined over Q
for G2 in ListG2 do
    C:=HyperellipticCurveFromG2Invariants(G2);
    if BaseField(C) eq RationalField() then ReducedWamelenModel(C); end if;
end for;
