//21/11/2018
// with the help of Xavier Xarles

jCM:={0,1728,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,-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};

load "FiberProductsDataadiconlymaximal.m";
k:=["4X7","5Nn"];
X:=FiberProducts[k]`curve;

printf "##########################\n %o \n##########################\n",k;
//------------------------------------------------
	ptsX:= Points(X: Bound:=10^6);

	G1:=k[1];
	G2:=k[2];
	j1:=Groups[G1]`jmap;	
	j2:=Groups[G2]`jmap;
	
	Pts:=[pt : pt in ptsX | (Evaluate(Denominator(j1),pt[1]) ne 0) and (Evaluate(Denominator(j2),pt[2]) ne 0) and (pt[3] ne 0)];
	assert &and{Evaluate(j1,pt[1]) eq  Evaluate(j2,pt[2]) : pt in Pts};	
	js:={Evaluate(j1,pt[1]): pt in Pts};
	for pt in Pts do
		printf "%o -> %o\n",pt,Evaluate(j1,pt[1]);
	end for;
    printf "-----------------------> %o\n", js diff jCM;

//------------------------------------------------
	
_,H1,mp1 := IsHyperelliptic(X);
C,mp2 := SimplifiedModel(H1);
F<x>,G<x> := HyperellipticPolynomials(C);
J := Jacobian(C);
rJ:=RankBound(J);	

boo,map:=IsIsomorphic(X,C);
assert boo;

printf " Genus = 2\n Hyperelliptic model: y^2=%o\n Rank Jacobian = %o\n",F,rJ;
 
ptsC := Points(C : Bound:=10^3);

printf " Points Search: \n %o\n",ptsC;

S:={pt[1]/pt[3] : pt in ptsC |pt[3] ne 0} join {pt[1] : pt in ptsC |pt[3] eq 0};
Hk,AtoHk:=TwoCoverDescent(C:PrimeBound:=30);
assert #S eq #Hk;
g:=(-2 + 3*x - 4*x^2 + 4*x^3);
assert IsDivisibleBy(F,g);
K<a>:=NumberField(g);
P<x> := PolynomialRing(K);
Fg:=Factorization(-2 + 3*x - 4*x^2 + 4*x^3);
h1:=Fg[1][1];h2:=Fg[2][1];
// h1; h2;
f:=h2*(x^2+1);

printf " Cover: \n y1^2=(x-2)*%o,\n y2^2=%o\n",  h1,f;

P1:=ProjectiveSpace(Rationals(),1);

S0:={pt[1]/pt[3] : pt in ptsC | pt[3] ne 0};

S:= {[d0,Evaluate(f, d0)] : d0 in S0 | not IsSquare(Evaluate(f, d0))};
Lpts:={};

printf " Twists {1} union %o \n",{s[2] : s in S};

// delta=1

	H:=HyperellipticCurve(f);
	HtoP1:=map<H->P1|[H.1,H.3]>;
	E,m:=EllipticCurve(H);
	EToP1:=Expand(Inverse(m)*HtoP1);
	MWgrp,MWmap,success,suc2:=MordellWeilGroup(E);
	//success; suc2;
	assert success;
	assert suc2;

printf " =====\n [oo] delta=1\n MW = %o\n ",	AbelianInvariants(MWgrp);

	V,R:=Chabauty(MWmap,EToP1);
	//V;
	//R;

	L1:={ EToP1( MWmap( v ) ) : v in V };

printf "Elliptic Chabauty: %o\n",L1;

	Lpts:=Lpts join L1;
	
// delta=/=1

for dd in S do
	xdelta:=dd[1];
	delta:=dd[2];
	Hdelta:=HyperellipticCurve(delta*f);	
	Pdelta:=Hdelta![xdelta,delta];	
	HdeltatoP1:=map<Hdelta->P1|[Hdelta.1,Hdelta.3]>;
	Edelta,m:=EllipticCurve(Hdelta,Pdelta);	
	EdeltaToP1:=Expand(Inverse(m)*HdeltatoP1);
	success,MWgrp,MWmap:=PseudoMordellWeilGroup(Edelta);
	//success; suc2;
	assert success;
	assert suc2;
printf " =====\n [%o] delta=%o -- \n MW = %o\n ",xdelta,delta, AbelianInvariants(MWgrp);
	V,R:=Chabauty(MWmap,EdeltaToP1:IndexBound:=2);
	Ldelta:={ EdeltaToP1( MWmap( v ) ) : v in V };
printf "Elliptic Chabauty: %o\n",Ldelta;
	Lpts:=Lpts join Ldelta;
end for;

	
	ptsC2:=&join{Points(C,xo[1]) : xo in Lpts | xo[2] ne 0} join PointsAtInfinity(C);;
	
	assert ptsC eq ptsC2;
	
	printf "\n\n SOLUTION:\n\n C(Q)=%o\n\n", ptsC2;


//------------------------------------------------
	ptsX2:={};
	for p in ptsC do
		ptsX2:={pt : pt in RationalPoints(p @@ map)} join ptsX2;
	end for;
	Pts2:=[pt : pt in ptsX2 | (Evaluate(Denominator(j1),pt[1]) ne 0) and (Evaluate(Denominator(j2),pt[2]) ne 0) and (pt[3] ne 0)];
	assert &and{Evaluate(j1,pt[1]) eq  Evaluate(j2,pt[2]) : pt in Pts2};	
	js:={Evaluate(j1,pt[1]): pt in Pts2};
	for pt in Pts2 do
		printf "%o -> %o\n",pt,Evaluate(j1,pt[1]);
	end for;

	assert ptsX eq ptsX2;

//  printf "-----------------------> %o\n", js diff jCM;
//------------------------------------------------
	printf "\n\n Rational points on the modular curve associated to %o:\n\n %o \n\n", k,ptsX2;
	printf "\n\n corresponding j-invariants:\n\n %o \n\n", js;
	printf "\n Non-CM j-invariants:\n %o\n", js diff jCM;


quit;
