SetAutoColumns(false);
SetColumns(100);
SetEchoInput(true);
Attach("routines.m");
Qx<x>:=PolynomialRing(Rationals());

P2<u,v,w>:=ProjectiveSpace(Rationals(),2);
Q1:=-4*u^2 - 4*v*w + 4*w^2;
Q2:=-2*u^2 + 2*u*w - 4*v*w + 2*w^2;
Q3:=2*u^2 + 4*u*v + 4*v^2;

Factorisation(BadPrimes([Q1,Q2,Q3]));
//Only bad primes are 2,5. That leaves 1,-1,2,-2,5,-5,10,-10 as possible
//values for delta.

//In fact, there are more bad primes for C than there are where Q1,Q2,Q3
//simultaneously vanish:
Factorisation(BadPrimes(DefiningEquations(SingularSubscheme(Curve(P2,Q1*Q3-Q2^2)))));

d:=-1;
pi,xpolmap,J:=Prym(P2,d*Q1,d*Q2,d*Q3);
D<U,V,W,R,S>:=Domain(pi);
Dpl:=Curve(Projection(Projection(D)));
//verify that negative d lead to D(R) empty.
assert HasRealPoints(Dpl) eq false;

d:=2;
pi,xpolmap,J:=Prym(P2,d*Q1,d*Q2,d*Q3);
D<U,V,W,R,S>:=Domain(pi);
assert not(IsLocallySolvable(D,2:AssumeIrreducible,AssumeNonsingular));

d:=5;
pi,xpolmap,J:=Prym(P2,d*Q1,d*Q2,d*Q3);
D<U,V,W,R,S>:=Domain(pi);
assert not(IsLocallySolvable(D,5:AssumeIrreducible,AssumeNonsingular));

d:=10;
pi,xpolmap,J:=Prym(P2,d*Q1,d*Q2,d*Q3);
D<U,V,W,R,S>:=Domain(pi);
assert not(IsLocallySolvable(D,2:AssumeIrreducible,AssumeNonsingular));

//So the only value left is delta=1:
d:=1;
pi,xpolmap,J:=Prym(P2,d*Q1,d*Q2,d*Q3);
D<U,V,W,R,S>:=Domain(pi);
P4:=Ambient(D);
C:=Codomain(pi);

//there are rational points on C and D:
p0:=C![0,1,0];
P0:=D!RationalPoints(p0@@pi)[1];P0;



KJ:=KummerSurface(J);
Km:=Scheme(KJ);
P3<k1,k2,k3,k4>:=Ambient(Km);

Ckm:=KummerImage(pi,xpolmap,J);

//verify that J(Q) is of rank 1 and has no torsion
assert TwoSelmerGroupData(J) eq 1 and TorsionBound(J,10) eq 1;

//Determine a generator of the Mordell-Weil group: First search for some free
//point
pt:=ReducedBasis(RationalPoints(J:Bound:=100))[1];
//Search up to the height bound suggested by that point:
RationalPoints(J:Bound:=Floor(CanonicalHeight(pt)+
                  Exp(HeightConstant(J:Effort:=100,Factor))));
//We see we can take:
g:=[elt<J|x^2 + 4*x - 4, 68*x - 64, 2>];

//As we'll see, a planar model of C will come in handy.
P2K<K1,K2,K3>:=ProjectiveSpace(BaseRing(Km),2);
Cpl:=P2K!!Projection(Ckm,P3![0,0,0,1]);
Cpleq:=DefiningEquation(Cpl);
Cpleq:=LCM([Denominator(c):c in Coefficients(Cpleq)])*Cpleq;
PolynomialRing(Integers(),3)!Cpleq;

//Create the Mordell-Weil group of J abstractly, together with its realisation
G:=AbelianGroup([0]);
mwmap:=map<G->J|a:->&+[Eltseq(a)[i]*g[i]: i in [1]]>;

p:=13;

//compute the reductions of everything mod 13
quot,emb:=Reduction(mwmap,p);
quot:=hom<Domain(quot)->Image(quot)|
     [quot(g):g in OrderedGenerators(Domain(quot))]>;
Kmp:=BaseChange(Km,Coercion(Rationals(),GF(p)));
Ckmp:=Scheme(Ambient(Kmp),
    DefiningEquations(BaseChange(Ckm,Coercion(Rationals(),GF(p)))));
Cplp:=BaseChange(Cpl,Coercion(Rationals(),GF(p)));
Jp:=Codomain(emb);
KJp:=KummerSurface(Jp);
mp:=func<v|
     Scheme(Ckmp,
        [v[1]*Ckmp.2-v[2]*Ckmp.1,
         v[2]*Ckmp.3-v[3]*Ckmp.2,
         v[3]*Ckmp.1-v[1]*Ckmp.3])>;
//We compute Ckmp(GF(13)) by taking Cplp(GF(13)) (the planar model!)
//and taking the GF(13) points above them on Ckmp.
V:=&join[RationalPoints(mp(v)): v in RationalPoints(Cplp)];
print "#of points on Ckmp:",#V;

//We lift the points to Jp(GF(13)) and keep one representative for each.
V:=&join[ #W eq 0 select {} else {W[1]} where
               W:=RationalPoints(Jp,KJp!Eltseq(v)): v in V];
print "#of points on Jacobian (divided by 2):",#V;

//pull points back to abstract group Jp(GF(13))
V:={v@@emb: v in V};
//take only those points in the image of the Mordell-Weil group
V:={Image(quot)!v: v in V | v in Image(quot)};
//include -points as well, since D is symmetric.
V:=V join {-v:v in V};
print "#of intersections with MWgroup:",#V;

//verify it's the generator and - the generator of the mordell-weil group
//that map to these points.
assert {quot(G.1),quot(-G.1)} eq V;
//verify that mod 13, the generator has indeed order 10.
assert Order(quot(G.1)) eq 10;
//Verify that the generator indeed maps onto Ckm:
Ckm!Eltseq(KJ!mwmap(G.1));

//Here, we are lazy: Of course, we should do these computations with limited
//11-adic information. However, numbers are still small enough in this case to do
//the computation over Z.
crd:=Eltseq(KJ!mwmap(11*G.1));

//Get the equations defining Ckm evaluated at crd.
L:=[Integers()|Evaluate(f,crd):f in DefiningEquations(Ckm)];

//indeed, mod 13^2, 11*G.1 does not map onto Ckm.
[L[1] mod 13^2, L[2] mod 13^2];
