
AttachSpec("echidna-2.0/src/echidna.spec");
load "LuerothInv.m";

//The goal of these algorithms is to generate quartics in one of the connected components  L1  in the intersection  L ^ D  .



function GenerateCubic(N0);
//Input : A randomness parameter  N0  .
//Output : A cubic surface  S  , its equation  Seq  , and two skew lines  l,m  on  S  .
//         (The lines  l,m  are both given by two of their points.)

stop := false;

while not stop do

//Ambient spaces:
P2<x,y,z> := ProjectiveSpace(Rationals(),2);
R2<x,y,z> := CoordinateRing(P2);

//Creating what will be the six random points:
Pts := [];

//Using the action of  SL_3  , we can take the first four points as follows:
Append(~Pts, [0,0,1]);
Append(~Pts, [0,1,0]);
Append(~Pts, [1,0,0]);
Append(~Pts, [1,1,1]);
//Random points:
a := Rationals() ! Random([-N0..N0]);
b := Rationals() ! Random([-N0..N0]);
c := Rationals() ! Random([-N0..N0]);
Append(~Pts, [a,b,c]);
a := Rationals() ! Random([-N0..N0]);
b := Rationals() ! Random([-N0..N0]);
c := Rationals() ! Random([-N0..N0]);
Append(~Pts, [a,b,c]);

//Determining the linear system:
Mons2 := MonomialsOfDegree(R2,3);
ME := [];
for pt in Pts do
    h := hom<R2 -> Rationals() | [pt[1],pt[2],pt[3]]>;
    for mon in Mons2 do
        Append(~ME,h(mon));
    end for;
end for;
M := Matrix(Rationals(),#Pts,#Mons2,ME);
N := NullSpace(Transpose(M));
d := Dimension(N);

//Non-degeneracy condition:
if d eq 4 then

//Creation of the Clebsch morphism:
cs := [];
for i:=1 to d do
    c := R2 ! 0;
    v := N.i;
    for i:=1 to #Mons2 do
        c := c + v[i]*Mons2[i];
    end for;
    Append(~cs,c);
end for;

//We have the morphism. Now to find its image.
P3<X,Y,Z,W> := ProjectiveSpace(Rationals(),3);
R3<X,Y,Z,W> := CoordinateRing(P3);
c := hom<R3 -> R2 | cs>;
Mons3 := MonomialsOfDegree(R3,3);
Mons3Ims := [];
for mon in Mons3 do
    Append(~Mons3Ims,c(mon));
end for;
Mons := MonomialsOfDegree(R2,9);
ME := [];
for monim in Mons3Ims do
    for mon in Mons do
        Append(~ME,MonomialCoefficient(monim,mon));
    end for;
end for;
M := Matrix(Rationals(),#Mons3Ims,#Mons,ME);
N := NullSpace(M);
v := N.1;
SEq := R3 ! 0;
for i:=1 to #Mons3 do
    SEq := SEq + v[i]*Mons3[i];
end for;

//Non-singularity condition for the cubic surface:
S := Scheme(P3,SEq);
d := Dimension(SingularSubscheme(S));
stop := (d eq -1);

end if;

end while;

//Creation of the skew lines:
p121 := [0,2/3,1/3];
p122 := [0,1/3,2/3];
p131 := [2/3,0,1/3];
p132 := [1/3,0,2/3];
P121 := [];
for c in cs do
    Append(~P121,Evaluate(c,p121));
end for;
P122 := [];
for c in cs do
    Append(~P122,Evaluate(c,p122));
end for;
P131 := [];
for c in cs do
    Append(~P131,Evaluate(c,p131));
end for;
P132 := [];
for c in cs do
    Append(~P132,Evaluate(c,p132));
end for;
p1 := P121;
p2 := P122;
q1 := P131;
q2 := P132;

return S,SEq,[p1,p2],[q1,q2];

end function;



function EquationsForLine(ps);
//Input : Two points in a projective space of dimension  3  .
//Output : Two equations determining line passing through these points in dual coordinates.

K := Parent(ps[1][1]);
M := Matrix(K,2,4,ps[1] cat ps[2]);
N := NullSpace(Transpose(M));
Es := [Eltseq(N.1), Eltseq(N.2)];

return Es;

end function;



function Degree2FromRamificationDivisor(P1,D);
//Input : A projective line  P1  and a degree  2  divisor  D  on  P1  .
//Output : A morphism over the base field of  P1  ramifying over  D  .

R<x,y> := CoordinateRing(P1);
S := Support(D);

//Distinguishing between rational and non-rational cases:
if #S eq 2 then
    B1 := Basis(Ideal(S[1]))[1];
    B2 := Basis(Ideal(S[2]))[1];
    f := map<P1 -> P1 | [B1^2,B2^2]>;
else
    B1 := Basis(Ideal(S[1]))[1];
    a := MonomialCoefficient(B1,x^2);
    b := MonomialCoefficient(B1,x*y);
    c := MonomialCoefficient(B1,y^2);
    t := b/a;
    n := c/a;
    //Distinguishing between the case with trivial trace and the generic case:
    if t eq 0 then
        f := map<P1 -> P1 | [x^2 - 2*n*x*y - n*y^2, x^2 + 2*n*x*y - n*y^2]>;
    else
        f := map<P1 -> P1 | [t*x^2 + 2*(t^2 - 2*n)*x*y + (t^3 - 3*n*t)*y^2, t*x^2 + 4*n*x*y + n*t*y^2]>;
    end if;
end if;

return f;

end function;



function CommonFiberBase(P1,f,fp);
//Input : A projective line  P1  and two degree  2  functions  f and fp  from  P1  to itself.
//Output : The unique element  q  of  P1  such that the fiber of  f  above  q  is also a fiber of  fp  .

R<x,y> := CoordinateRing(P1);
k := BaseRing(R);

//Recovering coefficients:
DEf := DefiningEquations(f);
DEfp := DefiningEquations(fp);
an := MonomialCoefficient(DEf[1],x^2);
bn := MonomialCoefficient(DEf[1],x*y);
cn := MonomialCoefficient(DEf[1],y^2);
ad := MonomialCoefficient(DEf[2],x^2);
bd := MonomialCoefficient(DEf[2],x*y);
cd := MonomialCoefficient(DEf[2],y^2);
apn := MonomialCoefficient(DEfp[1],x^2);
bpn := MonomialCoefficient(DEfp[1],x*y);
cpn := MonomialCoefficient(DEfp[1],y^2);
apd := MonomialCoefficient(DEfp[2],x^2);
bpd := MonomialCoefficient(DEfp[2],x*y);
cpd := MonomialCoefficient(DEfp[2],y^2);

//Some linear algebra:
M := Matrix(k,4,3,[an,bn,cn,-ad,-bd,-cd,apn,bpn,cpn,-apd,-bpd,-cpd]);
N := NullSpace(M);
v := N.1;

return [v[2],v[1]];

end function;



function InvolutoryPoint(SEq,ps,qs);
//Input : A pair of skew lines  l,m  on the cubic surface defined by  SEq  .
//Output : One of the involutory points  Q  for the pair  l,m  .

//Determining equations:
EsA := EquationsForLine(ps);
EsB := EquationsForLine(qs);

//Coordinatizing:
k := BaseRing(Parent(SEq));
R<X,Y,Z,W> := Parent(SEq);
P1<x,y> := ProjectiveSpace(k,1);
S<x,y> := CoordinateRing(P1);
hA := hom<R -> S | [x*ps[1][1] + y*ps[2][1], x*ps[1][2] + y*ps[2][2], x*ps[1][3] + y*ps[2][3], x*ps[1][4] + y*ps[2][4]] >;
hB := hom<R -> S | [x*qs[1][1] + y*qs[2][1], x*qs[1][2] + y*qs[2][2], x*qs[1][3] + y*qs[2][3], x*qs[1][4] + y*qs[2][4]] >;

//Map  p |--> T_p (S) ^ B  in coordinates:
M := Matrix(S,3,4,[hA(Derivative(SEq,X)),hA(Derivative(SEq,Y)),hA(Derivative(SEq,Z)),hA(Derivative(SEq,W)),
                    EsB[1][1],EsB[1][2],EsB[1][3],EsB[1][4],
                    EsB[2][1],EsB[2][2],EsB[2][3],EsB[2][4]]);
N := NullSpace(Transpose(M));
v := N.1;
M := Matrix(S,2,4,[qs[1][1],qs[1][2],qs[1][3],qs[1][4],
                   qs[2][1],qs[2][2],qs[2][3],qs[2][4]]);
v := Matrix(S,1,4,[S ! v[1], S ! v[2], S ! v[3], S ! v[4]]);
sol := Solution(M,v);
nf := sol[1][1];
df := sol[1][2];

//Map  q |--> T_q (S) ^ A  in coordinates:
M := Matrix(S,3,4,[hB(Derivative(SEq,X)),hB(Derivative(SEq,Y)),hB(Derivative(SEq,Z)),hB(Derivative(SEq,W)),
                    EsA[1][1],EsA[1][2],EsA[1][3],EsA[1][4],
                    EsA[2][1],EsA[2][2],EsA[2][3],EsA[2][4]]);
N := NullSpace(Transpose(M));
v := N.1;
M := Matrix(S,2,4,[ps[1][1],ps[1][2],ps[1][3],ps[1][4],
                   ps[2][1],ps[2][2],ps[2][3],ps[2][4]]);
v := Matrix(S,1,4,[S ! v[1], S ! v[2], S ! v[3], S ! v[4]]);
sol := Solution(M,v);
ng := sol[1][1];
dg := sol[1][2];

//Construction of the associated morphisms:
P1 := Curve(P1);
f := map<P1 -> P1 | [nf,df]>;
g := map<P1 -> P1 | [ng,dg]>;

//Branch divisor of  g  and a morphism ramifying over  it:
Rg := RamificationDivisor(g);
Bg := Pushforward(g,Rg);
fnew := Degree2FromRamificationDivisor(P1,Bg);

//Point below:
Q := CommonFiberBase(P1,f,fnew);
//Taking coordinates:
Q := [Q[1]*qs[1][1] + Q[2]*qs[2][1],Q[1]*qs[1][2] + Q[2]*qs[2][2],Q[1]*qs[1][3] + Q[2]*qs[2][3],Q[1]*qs[1][4] + Q[2]*qs[2][4]];

return Q;

end function;



function BranchQuartic(f,p);
//Input : A quaternary cubic form  f  and a point  p  on the corresponding surface
//Output : A quartic that is the ramification locus of the projection from  p  .

Rf<X,Y,Z,W> := Parent(f);
K := BaseRing(Rf);

//Normalizing the point  p  :
for i:=1 to 4 do
    if p[i] ne 0 then
       i0 := i;
       q := [p[1]/p[i],p[2]/p[i],p[3]/p[i],p[4]/p[i]];
    end if;
end for;

//Taking tangent directions in the induced coordinates:
R<x,y,z> := PolynomialRing(K,3);
S<t> := PolynomialRing(R);
teller := 1;
Mons := [x,y,z];
qvar := [];
for i:=1 to 4 do
    if i eq i0 then
        Append(~qvar, S ! q[i]);
    else
        Append(~qvar, S ! (q[i] + Mons[teller]*t));
        teller := teller + 1;
    end if;
end for;

//Determination of the quartic branch curve:
h := hom<Rf -> S | qvar>;
hf := h(f);
a := Coefficient(hf,3);
b := Coefficient(hf,2);
c := Coefficient(hf,1);
return b^2 - 4*a*c;

end function;



function GenerateL1(N0);
//Input : A randomness parameter  N0  .
//Output : A quartic  f  in  L1  .

//Putting everything together:
S,SEq,ps,qs := GenerateCubic(N0);
Q := InvolutoryPoint(SEq,ps,qs);
f := BranchQuartic(SEq,Q);
return f;

end function;












//Test:

Param := 2^4;
NumberOfQuartics := 2^6;
L1Quartics := [];

P2<x,y,z> := ProjectiveSpace(Rationals(),2);
R<x,y,z> := CoordinateRing(P2);
for i:=1 to NumberOfQuartics do
    i;
    f := GenerateL1(Param);
    Append(~L1Quartics,R!f);
end for;

for i:=1 to NumberOfQuartics do
    i;
    f := L1Quartics[i];
    I := DixmierOhnoInvariants(f);
    I3, I6, I9, J9, I12, J12, I15, J15, I18, J18, I21, J21, I27 := Explode(I);
    L := LuerothInv(I);
    if (I27 eq 0) and (L eq 0) then
        "Should be OK.";
    else
        "WARNING : Error in algorithm!";
    end if;
end for;

