/***
 *  Explicit computation for singular Luröth quartics
 *
 *  Distributed under the terms of the GNU Lesser General Public License (L-GPL)
 *                  http://www.gnu.org/licenses/
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation; either version 2.1 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *  Copyright 2012, R. Basson, R. Lercier, J. Sijsling & C. Ritzenthaler
 */

/***
 *
 * 
 * This small piece of MAGMA code aims at verifying that there is any 
 * degree 30 invariant vanishing on L2 and any degree 24 invariant 
 * vanishing on L1 (see [BLRS2012] for details). It is based on the
 * ECHIDNA package published by D. Kohel on his web page..
 * 
 * Bibliography
 * 
 * [BLRS2012] "An explicit expression of Lüroth invariant",  R. Basson,
 * R. Lercier, C. Ritzenthaler and J. Sijsling, 2012, preprint.
 *
 */ 

// Echidna package
AttachSpec("echidna-2.0/src/echidna.spec");


// Rational field
p := 0; FF:= Rationals();

// The algebra of invariants for quartics
W:= [3, 6, 9, 9, 12, 12, 15, 15, 18, 18, 21, 21, 27];
PJ<I3, I6, I9, J9, I12, J12, I15, J15, I18, J18, I21, J21, I27> := PolynomialRing(FF, W);



/* Random general quartics */
function SampleQuartic(: Pc := PolynomialRing(FF, 3), BD := 1)
    Base:= MonomialsOfWeightedDegree(Pc, 4);
    if Characteristic(FF) eq 0 then
	Coeff:= [ Random(-BD, BD) : i in [1..15] ];
    else
	Coeff:= [ Random(CoefficientRing(Pc)) : i in [1..15] ];
    end if;
    Coeff:= [ FF!Random(-BD, BD) : i in [1..15] ];

    return &+[Coeff[i]*Base[i] : i in [1..15]];
end function;



/* ____________________________________________________________________________*/
/* ____________________________________________________________________________*/


/* Degree 24 invariant for L1 ? */
"";"";"Degree 24 invariant for L1";"";

Margin := 50;
M24:= MonomialsOfWeightedDegree(PJ, 24);

/* Random Luröth quartics in L1 */    
load "L1Database.m";

/* Main computations
 ********************/

t := Cputime();

/* First matrix for general quartics */
MTRQ:= Matrix(FF,  #M24+Margin, #M24, []);
for i := 1 to #M24+Margin do
    form := SampleQuartic();
    J, W := DixmierOhnoInvariants(form);
    for j := 1 to #M24 do MTRQ[i][j] := Evaluate(M24[j], J); end for;
end for;

/* Second matrix for Luröth quartics in L1 */
MTRL:= Matrix(FF,  #M24+Margin, #M24, []);
for i := 1 to #M24+Margin do
    form := Random(D);
    J, W := DixmierOhnoInvariants(form);
    for j := 1 to #M24 do MTRL[i][j] := Evaluate(M24[j], J); end for;
end for;

"Matrices generated in", Cputime(), "s";
  
W := Vector(FF, [0 : i in [1..#M24+Margin]]);
sol, NsQ := Solution(Transpose(MTRQ), W);
Dimension(NsQ), "solution(s) found in", Cputime(t), "s";

W := Vector(FF, [0 : i in [1..#M24+Margin]]);
sol, NsL := Solution(Transpose(MTRL), W);
Dimension(NsL), "solution(s) found in", Cputime(t), "s";



/* ____________________________________________________________________________*/
/* ____________________________________________________________________________*/



/* Degree 30 invariant for L2 ? */
"";"";"Degree 30 invariant for L2";"";

Margin := 50;
M30:= MonomialsOfWeightedDegree(PJ, 30) diff {I3*I27};
/* Note that for a singular quartic I27 is zero ! */

/* Random Luröth quartics in L2 */
function SampleLurothL2(: Pc := PolynomialRing(FF, 3), BD := 4)    
    repeat 
	Coeff:= [ Random(-BD, BD) : i in [1..7] ];
	l1:= Pc.1 - Pc.2; l2:= Pc.2 - Pc.3; l3:= 2*Pc.1 - Pc.2 - Pc.3; l4:= Pc.1+Pc.2+Pc.3;
	l:= &+[ Coeff[i]*Pc.i : i in [1..3]];
	f:= l1*l2*l3*l4 + Coeff[4]*l1*l2*l3*l +  Coeff[5]*l1*l2*l4*l +
	    Coeff[6]*l1*l3*l4*l +  Coeff[7]*l2*l3*l4*l;
    until f ne 0;
    return f;
end function;



/* Main computations
 ********************/

t := Cputime();

/* First matrix for general quartics */
MTRQ:= Matrix(FF,  #M30+Margin, #M30, []);
for i := 1 to #M30+Margin do
    form := SampleQuartic();
    J, W := DixmierOhnoInvariants(form);
    for j := 1 to #M30 do MTRQ[i][j] := Evaluate(M30[j], J); end for;
end for;

/* Second matrix for Luröth quartics in L2 */
MTRL:= Matrix(FF,  #M30+Margin, #M30, []);
for i := 1 to #M30+Margin do
    form := Random(D);
    J, W := DixmierOhnoInvariants(form);
    for j := 1 to #M30 do MTRL[i][j] := Evaluate(M30[j], J); end for;
end for;

"Matrices generated in", Cputime(), "s";
  
W := Vector(FF, [0 : i in [1..#M30+Margin]]);
sol, NsQ := Solution(Transpose(MTRQ), W);
Dimension(NsQ), "solution(s) found in", Cputime(t), "s";

W := Vector(FF, [0 : i in [1..#M30+Margin]]);
sol, NsL := Solution(Transpose(MTRL), W);
Dimension(NsL), "solution(s) found in", Cputime(t), "s";
