// *****************************************************************************************************
// *** Relinearization Algorithm - Solve an overdetermined system with m >= (n+1)(n+2)/2-1=n*(n+3)/2 equations ***
// *** create a random polynomial quadratic system of m = (n+1)*(n+2)/2-1 equations in n variables; 
// *** the constant terms are chosen in such a way that the system is solvable
// *** the system is written out and read in again to be solved by the relinerazation method
// *** For your own work delete the first part and use your own file system.txt  
// *****************************************************************************************************


q:=7; // Parameters
GF:=GaloisField(q);
n:=3;
m:=Integers()!((n+1)*(n+2)/2-1);
Pol<[x]>:=PolynomialRing(GF,n);

pol:=[]; // create random system pol
for i:=1 to m do
  pol[i]:=Pol!0;
  for j:=1 to n do
    for k:=j to n do
      pol[i]:=pol[i]+Random(GF)*x[j]*x[k];
    end for;
    pol[i]:=pol[i]+Random(GF)*x[j];
  end for;
end for;


a:=Random(VectorSpace(GF,n));
pol1:=pol;  // ensure that the overdetermined system pol has a solution
for j:=1 to n do
  for i:=1 to m do
    pol1[i]:=Evaluate(pol1[i],x[j],a[j]);
  end for;
end for;


for i:=1 to m do
  pol[i]:=pol[i]-pol1[i];
end for;

printf"Write system.txt \n \n";
SetOutputFile("system.txt":Overwrite:=true);
printf "q:= %o ; \n \n", q;
printf "n:= %o ; \n \n", n;
printf "m:= %o ; \n \n", m;
printf "GF<w>:=GaloisField(q); \n \n";
printf "Pol<[x]>:=PolynomialRing(GF,n); \n \n";
printf "pol:= %o ; \n \n", pol;
UnsetOutputFile();

																

// **************************************************************************\
// solves the system pol using the Relinearization technique         ********\
// **************************************************************************

clear;
load "system.txt";

// compute the Macaulay matrix A of the system pol (without constant terms)
// compute the vector v containing the constant terms

A:=ZeroMatrix(GF,m,Integers()!((n+1)*(n+2)/2-1)); 

v:=VectorSpace(GF,m)!0;

for i:=1 to m do

  counter:=1;
  for j:=1 to n do // quadratic terms
    for k:=j to n do
	A[i][counter]:=MonomialCoefficient(pol[i],x[j]*x[k]);
	counter:=counter+1;
    end for;
  end for;
  for j:=1 to n do // linear terms
    A[i][counter]:=MonomialCoefficient(pol[i],x[j]);
    counter:=counter+1;
  end for;
  v[i]:= MonomialCoefficient(pol[i],1);     //constant terms;
end for;
 printf "Linear System: \n%o * sol = %o \n \n", A ,v;
 _,sollong:=IsConsistent(Transpose(A),(-v)); // solution of the linear system

sol:= Eltseq(sollong); // Take only the last n components of the solution \
				// (corresponding to the variables x1, to , xn
repeat
  sol:=Remove(sol,1);
until #(sol) eq n;

printf "solution:= %o \n \n", sol;

printf "Answer  see Relinearizationoutput.txt \n \n";

SetOutputFile("Relinearizationoutput.txt":Overwrite:=true);

printf "q:= %o ; \n \n", q;
printf "n:= %o ; \n \n", n;
printf "m:= %o ; \n \n", m;
printf "GF<w>:=GaloisField(q); \n \n";
printf "Pol<[x]>:=PolynomialRing(GF,n); \n \n";
printf "pol:= %o ; \n\n", pol;
printf "Linear System: \n%o * sol = %o \n\n", A ,v;
printf "solution:= %o \n\n", sol;

printf  "verify that we have found a valid solution\n\n";

for i:=1 to m do
  ans:= pol[i] ;
  for j:=1 to n do 
     ans := Evaluate(ans,x[j],sol[j]) ;
  end for ;
  printf "ans[%o]=%o\n",i,ans ;
end for ;

UnsetOutputFile();