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

P2<u,v,w>:=ProjectiveSpace(Rationals(),2);
Q1:=v^2+v*w-w^2;
Q2:=u^2-v^2-w^2;
Q3:=w^2+u*v;
Factorisation(BadPrimes([Q1,Q2,Q3]));
Factorisation(BadPrimes(DefiningEquations(SingularSubscheme(Curve(P2,Q1*Q3-Q2^2)))));

//We only have to consider delta in {1,-1,2,-2}. We do away with all of
//them except 1

d:=-1;
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:=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:=-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));

//that leaves one twist:

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);
KJ:=KummerSurface(J);
Km:=Scheme(KJ);
P3<k1,k2,k3,k4>:=Ambient(Km);

//We claim that D has points everywhere locally.
//First, it has real points:
Dpl:=Curve(Projection(Projection(D)));
Genus(Dpl);
HasRealPoints(Dpl);

for p in [p:p in [2..11]|IsPrime(p)] do
  if IsLocallySolvable(D,p:AssumeIrreducible,AssumeNonsingular) then
    print "Curve is solvable at prime ",p;
  end if;
end for;

for p in [p:p in [12..98] cat [656023]|IsPrime(p)] do
  Fp:=GF(p);
  Dp:=BaseChange(Dpl,Coercion(Rationals(),Fp));
  P2p:=AmbientSpace(Dp);
  repeat
    W:=RationalPoints(Dp meet Scheme(P2p,&+[Random(Fp)*P2p.i:i in [1..3]]));
    W:={P: P in W| IsNonSingular(Dp!P)};
  until #W gt 0;
  print "Curve is solvable at prime ",p;
end for;

//We determine the Mordell-Weil group of the Prym-variety J.
//The rank has an upper bound of 2
assert TwoSelmerGroupData(J) eq 2;
B:=ReducedBasis(RationalPoints(J:Bound:=1000));B;

//And to verify that this generates the entire Mordell-Weil group, notice that
//the height bound up to which we must search is already lower than
//what we have searched for. See: 
//M. Stoll, On the height constant for curves of genus two. II. Acta Arith. 104 ,
//2002, no. 2, 165--182; MR 2003f:11093.

Exp(HeightConstant(J:Effort:=1000,Factor)+
   CanonicalHeight(B[1])*CanonicalHeight(B[2])*
   Minimum(CanonicalHeight(B[1]+B[2]),CanonicalHeight(B[1]-B[2]))/(4*Re\
gulator(B)));

//compute the embedding of C in the Kummer surface
 
Ckm:=KummerImage(pi,xpolmap,J);

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

p:=11;
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:=BaseChange(Ckm,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(p)) by taking Cplp(GF(p)) (the planar model!)
//and taking the GF(p) points above them on Ckmp.
V:=&join[RationalPoints(mp(v)): v in RationalPoints(Cplp)];
V:=&join[ #W eq 0 select {} else {W[1]} where
             W:=RationalPoints(Jp,KJp!Eltseq(v)): v in V];
V:={v@@emb: v in V};
V:={Image(quot)!v: v in V | v in Image(quot)};
V:=V join {-v:v in V};

V11:=<quot,V>;

//this results in:
V11;

//lift this information to 11^2
psi:=DefiningEquations(Ckm)[2];
V11:=NextStage(V11,mwmap,psi,11,2);

V11;

//now do the same thing for 7

p:=7;
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:=BaseChange(Ckm,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])>;
V:=&join[RationalPoints(mp(v)): v in RationalPoints(Cplp)];
V:=&join[ #W eq 0 select {} else {W[1]} where
             W:=RationalPoints(Jp,KJp!Eltseq(v)): v in V];
V:={v@@emb: v in V};
V:={Image(quot)!v: v in V | v in Image(quot)};
V:=V join {-v:v in V};

V7:=<quot,V>;
V7;

CosetIntersection(V11,V7);
