SetAutoColumns(false);
SetColumns(100);
SetEchoInput(true);
Attach("routines.m");

P<x>:=PolynomialRing(Rationals());
f:=x*(x-1)*(x-2)*(x+1)*(x+3);
F:=HyperellipticCurve(f);
J:=Jacobian(F);
Km:=KummerSurface(J);
KJ:=Scheme(Km);
R<k1,k2,k3,k4>:=Parent(DefiningPolynomial(Km));

A<u,v,w>:=PolynomialRing(Rationals(),3);
L:=[];

//We compute the plane section k1+k2+t*k3+k4 for some
//rational values of t and determine the corresponding Q1,Q2,Q3
//from that. We store them with the corresponding value of t
for t in [t: t in [-20..20]| t notin [-2,0,2,3]] do
  print "t=",t;
  V:=k1+k2+t*k3+k4;
  
  //this routine does the hard work: It computes a model of the form
  //Q1*Q3=Q2^2 for V intersect Km.
  
  mp,Q1,Q2,Q3:=KummerSection(Km,V);
  c:=MonomialCoefficient(A!Q1,u^2);
  L cat:=[<t,A!Q1/c,A!Q2/c,A!Q3/c>];
end for;

//Now we pass to the function field.

Qt<t>:=RationalFunctionField(Rationals());
P<X>:=PolynomialRing(Qt);
f:=X*(X-1)*(X-2)*(X+1)*(X+3);
F:=HyperellipticCurve(f);
J:=Jacobian(F);
Km:=KummerSurface(J);
KJ:=Scheme(Km);
R<k1,k2,k3,k4>:=Parent(DefiningPolynomial(Km));
V:=k1+k2+t*k3+k4;
P2t<U,V,W>:=ProjectiveSpace(Qt,2);
mp:=hom<A->CoordinateRing(P2t)|[U,V,W]>;

// A routine that interpolates a rational function with numerator and
// denominator or degree at most d through the points described by L
// (your should really apply this to an overdetermined system)

function ratint(L,d)
  M:=[];
  for c in L do
    t0:=c[1];
    C:=c[2];
    M cat:=[[t0^i:i in [0..d]] cat [-C*t0^i:i in [0..d]]];
  end for;
  V:=Kernel(Transpose(Matrix(M)));
  if Dimension(V) gt 0 then
    print Basis(V);
    v:=Eltseq(Basis(V)[1]);
    return (Parent(x)!v[1..d+1])/(Parent(x)!v[d+2..2*d+2]);
  else
    error "no solution";
  end if;
end function;

//We now solve Q1,Q2,Q3 over Qt via interpolation:
mons:=MonomialsOfDegree(Parent(u),2);
cf:=[<l[1],[MonomialCoefficient(l[2],m):m in mons]>:l in L];
Q1t:=&+[ratint([<c[1],c[2][i]>:c in cf],5)*mp(mons[i]):i in [1..#mons]];

cf:=[<l[1],[MonomialCoefficient(l[3],m):m in mons]>:l in L];
Q2t:=&+[ratint([<c[1],c[2][i]>:c in cf],3)*mp(mons[i]):i in [1..#mons]];

cf:=[<l[1],[MonomialCoefficient(l[4],m):m in mons]>:l in L];
Q3t:=&+[ratint([<c[1],c[2][i]>:c in cf],3)*mp(mons[i]):i in [1..#mons]];

//a help function that computes the splitting discriminant of a singular
//plane conic.

function dsc(M)
  B:=Basis(Kernel(M));
  B:=ExtendBasis(B,Generic(Universe(B)));
  B:=Matrix([B[2],B[3]]);
  return -Determinant(B*M*Transpose(B));
end function;

//make the answers a bit nicer:

Q1:=36*(t-41/18)*Q1t;Q1;
Q2:=36*(t-41/18)*Q2t;Q2;
Q3:=36*(t-41/18)*Q3t;Q3;

//The curve is:
Ct:=Curve(P2t,Q1*Q3-Q2^2);Ct;

//compute the nodes of the Kummer surface:
nodes:=RationalPoints(SingularSubscheme(KJ));

//Here, we're lazy: We compute the tropes by taking all possible planes through
//3 nodes and only takiong those that have 6 nodes on them.

triples:=Subsets({a:a in nodes},3);
function plane(c)
  P:=Ambient(Scheme(Rep(c)));
  cf:=Basis(Kernel(Transpose(Matrix([Eltseq(p): p in c]))));
  return Scheme(P,[&+[v[i]*P.i:i in [1..Dimension(P)+1]]: v in cf]);
end function;
tropes:={v: c in triples | #{a: a in nodes | a in v} eq 6 where v:=plane(c)};

//Check that we have them all.
assert #tropes eq 16;

//Now compute the bitangents that these lines induce on Ct:
tng1t:=[Curve(P2t,Evaluate(DefiningEquation(trope),[U,V,W,-U-V-t*W])):trope in tropes];
[DefiningPolynomial(t):t in tng1t];

//These are the singular quadrics
sngquads:=[Q1-2*x*Q2+x^2*Q3: x in [0,-3,-1,1,2]] cat [Q3];
[Determinant(SymmetricMatrix(Q)):Q in sngquads];

//get the splitting discriminants:
dscs:=[dsc(SymmetricMatrix(Q)):Q in sngquads];dscs;

//Check that all split when one splits;
[IsSquare(l*dscs[1]):l in dscs];

//Therefore, it is sufficient to get the following to be a square:
dscs[1];

//parametrize the corresponding conic:
L:=DefiningEquations(Parametrization(Conic(HyperellipticCurve(
                            Numerator(dscs[1])))));

//and extract the desired substitution from that
sval:=Evaluate(L[1],[t,1])/Evaluate(L[3],[t,1]);sval;

//Create the induced maps:
scov:=hom<Qt->Qt| sval>;
lft:=hom<Parent(U)->Parent(U)| scov*Bang(Qt,Parent(U)), [U,V,W]>;

//Make sure that things are really printed with s now:
//(this ruins printing with t)
_<s>:=Codomain(scov);

tng1s:=[Scheme(P2t,lft(DefiningEquation(l))):l in tng1t];

//Take the primary components of the singular conics. These are bitangents as
//well.
tng2s:=&cat[PrimaryComponents(Scheme(P2t,lft(Q))):Q in sngquads];
Cs:= Curve(P2t,lft(DefiningEquation(Ct)));
tngs:=tng1s cat tng2s;

//we now really have 28 bitangents
#tngs eq 28;

//Make sure that things are really printed with s now:
_<s>:=Codomain(scov);

//Check that all these really are bitangents
[SingularSubscheme(I) eq ReducedSubscheme(I) where I:= Cs meet T: T in tngs];

//The defining equations of all 28 bitangents:
[DefiningEquation(T):T in tngs];
