/***
 *  Explicit computation of Luroth's invariant
 *
 *  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 computing Luroth's invariant,
 * over the rationals (see [BLRS2012] for details). Is 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, J. Sijsling and C. Ritzenthaler, 2012, preprint.
 *
 */ 

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

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

// Uncomment what follows for a finite field computation
p:= NextPrime(10^6); n:= 1; while p^n lt 2000 do n+:=1; end while; FF:= GF(p^n);

// 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);

M54:= MonomialsOfWeightedDegree(PJ, 54);

// Hadmard bound of a matrix
function HadamardBoundSquare(MTR)
    Nrm := 1; for i := 1 to Nrows(MTR) do
	Sm := 0;
	for j := 1 to Ncols(MTR) do
	    Sm +:= MTR[i,j]^2;
	end for;
	if Sm ne 0 then
	    Nrm *:= Sm;
	end if;
    end for;
    return Ceiling(Log(2, IntegerRing()!Nrm));
end function;
    
/* Random Luroth quartics */
function SampleLuroth1(: Pc := PolynomialRing(FF, 3), BD := 4)    
    repeat 
	if Characteristic(FF) eq 0 then
	    Coeff:= [ Random(-BD, BD) : i in [1..7] ];
	else
	    Coeff:= [ Random(CoefficientRing(Pc)) : i in [1..7] ];
	end if;
	Coeff:= [ FF!Random(-BD, BD) : i in [1..7] ];
	l1:= Pc.1; l2:= Pc.2; l3:= 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;

/* 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;

/* Dumping a matrix on a file */
procedure DumpMatrix(output, MTR)

    OutputFile := POpen("gzip > " cat output, "w");

    Puts(OutputFile, Sprintf("%o", Nrows(MTR)));
    Puts(OutputFile, Sprintf("%o", Ncols(MTR)));

    for i := 1 to Nrows(MTR) do
	for j := 1 to Ncols(MTR) do
	    Puts(OutputFile, Sprintf("%o", MTR[i,j]));
	end for;
    end for;

end procedure;

/* Restoring a matrix from a file */
function ReadMatrix(input : FF := IntegerRing())

    InputFile := POpen("gzip -d < " cat input, "r");

    str := Gets(InputFile); nr := StringToInteger(str);
    str := Gets(InputFile); nc := StringToInteger(str);

    MTR := Matrix(FF, nr, nc, []);

    for i := 1 to nr do
	for j := 1 to nc do
	    str := Gets(InputFile);
	    MTR[i,j] := StringToInteger(str);
	end for;
    end for;
    
    return MTR;

end function;

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

Margin := 50;
t := Cputime();

/* First matrix */
MTRQ:= Matrix(FF,  #M54+Margin, #M54, []);
for i := 1 to #M54+Margin do
    repeat
	form := SampleQuartic();
	J, W := DixmierOhnoInvariants(form);
    until &*J ne 0;
    for j := 1 to #M54 do MTRQ[i][j] := Evaluate(M54[j], J); end for;
end for;
"Estimated accurancy needed :", HadamardBoundSquare(MTRQ), "bits";  
DumpMatrix("MTRQ" cat IntegerToString(p) cat ".dat.gz", MTRQ);

/* Second matrix */
MTRL:= Matrix(FF,  #M54+Margin, #M54, []);
for i := 1 to #M54+Margin do
    repeat
	form := SampleLuroth1();
	J, W := DixmierOhnoInvariants(form);
    until &*J ne 0;
    for j := 1 to #M54 do MTRL[i][j] := Evaluate(M54[j], J); end for;
end for;
"Estimated accurancy needed :", HadamardBoundSquare(MTRL), "bits";     
DumpMatrix("MTRL" cat IntegerToString(p) cat ".dat.gz", MTRL);

"Matrices generated in", Cputime(), "s";

/* Kernel of the first matrix */
W := Vector(FF, [0 : i in [1..#M54+Margin]]);
sol, NsQ := Solution(Transpose(MTRQ), W);
Dimension(NsQ), "solution(s) found in", Cputime(t), "s";
if Dimension(NsQ) ne 215 then error "Mauvaise dimension"; end if;
DumpMatrix("NsQ" cat IntegerToString(p) cat ".dat.gz", Matrix(Basis(NsQ)));
//sub<VectorSpace(FF, Ncols(MTRL))|[Matrix(Basis(NsQ))[i] : i in [1..Nrows(Matrix(Basis(NsQ)))]]> eq NsQ;

/* Kernel of the second matrix */
W := Vector(FF, [0 : i in [1..#M54+Margin]]);
sol, NsL := Solution(Transpose(MTRL), W);
Dimension(NsL), "solution(s) found in", Cputime(t), "s";
if Dimension(NsL) ne 216 then error "Mauvaise dimension"; end if;
DumpMatrix("NsL" cat IntegerToString(p) cat ".dat.gz", Matrix(Basis(NsL)));

/* Luroth invariant */
L:= Eltseq(Basis(Complement(NsL, NsQ))[1]);
Sol:= &+[ M54[i]*L[i] : i in [1..#M54] ];
Write("LurothChar" cat IntegerToString(p) cat ".m", "L :=" : Overwrite:= true);
Write("LurothChar" cat IntegerToString(p) cat ".m", Sol : Overwrite:= false);
Write("LurothChar" cat IntegerToString(p) cat ".m", ";" : Overwrite:= false);

/* Result */
Sol;
