// load "sparse_issac.magma";

// different benchmarks [n,t,nb] where n is the number of unknowns, t the size or the type of the support and nb the number of equations


// Small example: Fewnomials very sparse
EX1:=[60,60+3,59];

// Medium size example: Fewnomials sparse
EX2:=[150,3*150,434];

// Bilinear system: two blocks of variables of size (7,30)
EX3:=[37,-7,76];

// Choose one benchmark
bench:=EX2;
n:=bench[1];t:=bench[2];nb:=bench[3];
force:=1; // force the random system to have one solution (optional)
p0:=65521; // size of the finite field

// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------

K:=GF(p0);
R<[x]>:=PolynomialRing(K,n,"grevlex");
//assigns nice names to the variables
if t lt 0 then
    AssignNames(~R,["x" cat IntegerToString(-t-i): i in [1..(-t)]] cat ["y" cat IntegerToString(n+t-i) : i in [1..(n+t)]]);
    H1:=x[-t];H2:=x[n];
else
    AssignNames(~R,["x" cat IntegerToString(i): i in [1..(n-1)]] cat ["h"]);
    H1:=x[n];H2:=x[n];
end if;

// support of polynomial system is of size t + admissible signatures
load "data.magma";

// Random coefficients
eqs:=[&+[Random(K)*u : u in support] : i in [1..nb]];

// this is to ensure that we have a solution (optional)
if force eq 1 then
    sol:=[Random(K) : i in [1..(n-1)]] cat [1];
    eqs:=[h-Evaluate(h,sol)*H1*H2/Evaluate(H1*H2,sol) : h in eqs];
end if;

// There is a bug in magma ...
printf "sparse gb in degree 1 ";
time gb0:= GroebnerBasis(GroebnerBasis(eqs,2),2);

if #gb0 ne nb then
    printf "****************************** Equations are not linearly independent !!\n";
end if;

// list of admissible signatures:  if Fsign[i] = [j1, j2 , ...] then support[j1]*F[i],support[j2]*F[i], ... are admissible
Fsign:=[u cat sign : u in Fsign];

// polynomial representation of the matrix
F:=&cat[[support[j]*gb0[i] : j in Fsign[i]] : i in [1..nb]];

// linear algebra part: first we build the matrix
printf "Build A\n";

mi:=&join [{u : u in Monomials(h)} : h in F];
SortMon:=function(a,b) if a lt b then return 1; else return -1; end if; end function;
mi:=Sort([u : u in mi],SortMon);
N:=#F; // number of rows
M:=#mi; // number of columns

Mons:=AssociativeArray();
for k:=1 to M  do
    Mons[mi[k]]:=k;
end for;

A:=Matrix(K,N,M,[0 : i in [1..N*M]]);
for i in [1..N] do
    for m in Terms(F[i]) do
	cm:=LeadingCoefficient(m);
	A[i,Mons[m/cm]]:=cm;
    end for;
end for;
printf "done\n";

// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------

printf "sparse gb in degree 2 (Echelon Form of a matrix of size %ox%o) ",N,M;
time B:=EchelonForm(A);
// display the 10 last elements
FF:=&cat[ [&+[B[i,j]*mi[j]: j in [1..M]]] : i in [(M-10)..(M-1)]];

printf "Non Reduced Grobner Basis (last 10 elements) %o .... \n",FF;

printf "Affine Solutions (last 10 elements) %o .... \n",[NormalForm(h,[H2-1,H1-1]) : h in FF];

//----------------------------------------------------------------------
printf "and now compute the GB with magma --> ";
time gb:=GroebnerBasis(eqs);