// *****************************************************************************************************
// *** create a random polynomial quadratic system of m  equations in n variables; 
// *** for solving the system by the method of section 8.7.3  one needs  nu=n/m >= 2
// *** the system is written out and read in again 
// *** For your own work delete the first part and use your own file system.txt  
// **************************************************************************************************** 
clear() ;

q:=4;
GF<w>:=GaloisField(q);
m:=3;
n:=7;

printf "******************************************************************** \n";
printf "*** Generate Random System of m equations in n=nu*m variables    *** \n";
printf "*** nu>=2 is required                                            *** \n";
printf "*** given is nu=%o                                            ***\n",n/m;
printf "******************************************************************** \n\n\n";

printf "*** Parameters ***\n\n";
printf "GF:=GaloisField(%o) \n\n", q;
printf "m:= %o \n\n", m;
printf "n:= %o \n\n\n", n;
printf "nu=%o\n\n",n/m;

Pol<[x]>:=PolynomialRing(GF,n);

pol:=[];    // random quadratic polynomial
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;
    pol[i]:=pol[i]+Random(GF);
end for;

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



load "system.txt";

PoltoMat:= function(m,n,pol) // transforms the system pol into matrix form
  Mat:=[];
  lin:=[] ;
  con:=[];
  for i:=1 to m do
    Mat[i]:=ZeroMatrix(GF,n,n);
    lin[i]:=ZeroMatrix(GF,n,1) ;
    con[i]:=MonomialCoefficient(pol[i],1);
    for j:=1 to n do
      lin[i][j][1]:=MonomialCoefficient(pol[i],x[j]) ;
      for k:=j to n do
         Mat[i][j][k]:=MonomialCoefficient(pol[i],x[j]*x[k]);
        end for;
    end for;
  end for;
  return Mat,lin,con;
end function;


MattoPol:= function(mat,lin,con,m,n) // transforms the matrix form into polynomials
  poll := [] ;
  for i:=1 to m do
    poll[i] := Pol!con[i];
    for j:=1 to n do
      poll[i] +:= lin[i][j][1]*x[j] ;
      for k:=1 to n do
        poll[i] +:=  mat[i][j][k]*x[j]*x[k];
      end for ;
    end for;
  end for;
  return poll;
end function;

Mat,Lin,Con:=PoltoMat(m,n,pol);  // compute the matrices representing the homogeneous quadratic part of the system pol 

orig:=pol ;

//SetOutputFile("answer.txt":Overwrite:=true);

nu := n/m ;
omega := Floor(nu) ;

printf "nu = %o   omega = %o\n\n",nu,omega ;
printf "pol:= \n%o \n\n", pol;
printf "Mat:= \n%o \n\n", Mat;
 
T:=[];
Ta:=[];

 // Algorithm 8.11
printf "Find a linear transformation transforming the system pol into the form of Figure 8.2 \n";

for loop:=2 to m do  // main loop
printf " --------------- %o th step ----------------------------- \n\n", loop-1;
  k:=Minimum(omega-1, loop-1);
  printf"k:=%o \n \n",k;
  V:=VectorSpace( GF,k*m);

// Step 3
  T[loop]:=IdentityMatrix(Pol,n); // T[loop] is a matrix containing unknowns in the loop-th row
  
  for i:=1 to n do   
    T[loop][i][loop]:=x[i];
  end for;

  T[loop][loop][loop]:=1;

  printf"T[%o]= \n%o\n\n",loop, T[loop];
  
  A:=[];
  v:=[];
   
  for i:=1 to m do  // Compute the unknowns in T[loop] in such a way that
    M:= Transpose(T[loop]) * Mat[i] * T[loop];   
    for j:=1 to k do
      a:= M[j][loop]+M[loop][j];
      for t:=1 to n do
        A := A cat [MonomialCoefficient(a,x[t])];
      end for;
      v:=v cat [MonomialCoefficient(a,1)];
    end for;
  end for;
    
  AMat:=Matrix( GF,k*m,n,A); 
  v:=V!v;

// Step 4
  _,vec,ker:=IsConsistent(Transpose(AMat),v);
//  vec:=vec+Random(ker); //  

  printf "A= \n%o \n \n", AMat;
  printf "v= %o \n \n", v;
  printf "vec= %o \n \n", vec;

// Step 5 
  Ta[loop]:=IdentityMatrix( GF,n);  //Transformation Matrix
  for i:=1 to n do
    Ta[loop][i][loop]:=vec[i];
  end for;
  Ta[loop][loop][loop]:=1;

  for i:=1 to m do
    Mat[i]:=Transpose(Ta[loop]) * Mat[i] * Ta[loop];   //Transform the system
  end for;

end for; // end of main loop


TMat:=IdentityMatrix( GF,n); //Step 7
for i:=2 to m do
  TMat:=TMat*Ta[i];
end for;

// Compute simplified system ---------------------------

Tpol:=[];   // transform TMat into a polynomial (makes it easier to deal with the linear part of the system)
for i:=1 to n do
    Tpol[i]:=Pol!0;
    for j:=1 to n do
        Tpol[i]:=Tpol[i]+TMat[i][j]*x[j];
    end for;
end for;

// Linear Transformation of the original system

Q:=pol;
for loop:=1 to m do
    for i:=1 to n do
        Q[loop]:=Evaluate(Q[loop],x[i],Tpol[i]);
    end for;
end for;

printf "Simplified System= \n %o \n", Q;    // This is the system of Figure 8.2


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

Qorg:=Q;
// Now starts Algorithm 8.12

QMat:=PoltoMat(m,n,Q);     // Matrices representing the homogeneous quadratic part
printf "QMat:= %o \n\n", QMat;

//Step 1: Set the linear equations L_1, \dots , L_{omega-1} (x_m+1, \dots, x_n) to zero
A:=ZeroMatrix( GF,(omega-1)*m,n-m);
for i:=1 to m do
  for k:=1 to omega-1 do
    for j:=1 to n-m do
      A[(omega-1)*(i-1)+k][j]:=MonomialCoefficient(Q[i],x[k]*x[j+m]);
    end for;
  end for;
end for;

v:=[];
for i:=1 to m do
  for k:=1 to omega-1 do
    v[(omega-1)*(i-1)+k]:=MonomialCoefficient(Q[i],x[k]);
  end for;
end for;
Vm:=VectorSpace( GF,m*(omega-1));
v:=Vm!v;

printf"A=\n %o \n\n",A;
printf"v= %o \n\n",v;

   printf "Rank(A)=%o full rank would be %o\n\n",Rank(A),m;

_,sollin,ker:= IsConsistent(Transpose(A),v);
sollin:=sollin+Random(ker);
printf "sollin= %o \n \n", sollin; // solution of the linear equations = values of the variables x_{m+1}, \dots, x_n

for i:=1 to m do  // Step 2: substitute the values of x_{m+1}, \dots, x_n into the system
  for j:=1 to n-m do
    Q[i]:=Evaluate(Q[i],x[j+m], sollin[j]);
  end for;
end for;
Qdet:=Q;
printf "Simplified System= \n %o \n", Q;  // System in x_1, \dots, x_m

//Step 3
M:=ZeroMatrix( GF,m,omega-1+Integers()!((m-omega+2)*(m-omega+3)/2));
for i:=1 to m do
  for j:=1 to omega-1 do
    M[i][j]:=MonomialCoefficient(Q[i],x[j]^2);
  end for;
  counter:=omega;
  for j:=omega to m do
    for k:=j to m do
      M[i][counter]:=MonomialCoefficient(Q[i],x[j]*x[k]);
      counter:=counter+1;
    end for;
    M[i][counter]:=MonomialCoefficient(Q[i],x[j]);
    counter:=counter+1;
  end for;
  M[i][counter]:=MonomialCoefficient(Q[i],1);
end for;

printf "M: \n%o \n\n", M;
Mech, Sech:=EchelonForm(M);
printf "Mech: \n%o \n\n", Mech;

                            
Qech:=[];  // Echelon Form of the determined system
for i:=1 to m do;
  Qech[i]:=Pol!0;
  for j:=1 to omega-1 do
    Qech[i]:=Qech[i]+Mech[i][j]*x[j]^2;
  end for;
  counter:=omega;
  for j:=omega to m do
    for k:=j to m do
      Qech[i]:=Qech[i]+Mech[i][counter]*x[j]*x[k];
      counter:=counter+1;
    end for;
    Qech[i]:=Qech[i]+Mech[i][counter]*x[j];
    counter:=counter+1;
  end for;
  Qech[i]:=Qech[i]+Mech[i][counter];
end for;
Qechorg:=Qech;

printf"Qech:\n%o \n\n", Qech;

a:=omega-1;

Poly<[y]>:=PolynomialRing( GF,m-omega+1,"grevlex");    // Polynomial Ring in m-omega+1 variables in which the Groebner basis is computed

pol2:=[];      // build system
for loop:=1 to m-omega+1 do
  pol2[loop]:=Poly!0;
  for i:=1 to m-omega+1 do
    for j:=i to m-omega+1 do
      pol2[loop]:=pol2[loop]+MonomialCoefficient(Qech[a+loop],x[a+i]*x[j+a])*y[i]*y[j];
    end for;
    pol2[loop]:=pol2[loop]+MonomialCoefficient(Qech[loop+a],x[i+a])*y[i];
  end for;
  pol2[loop]:=pol2[loop]+MonomialCoefficient(Qech[loop+a],1);
end for;

printf"pol2= %o \n\n", pol2; ////////////////////////////////////////////////////////////////////////

I:=Ideal(pol2); // Step 5: Solve the quadratic system
V:=Variety(I); 
V;             // Solution in x_{m-omega, \dots, x_m}

for loop:=1 to m do  // Substitute the solution into Qech
  for i:= omega to m do
    Qech[loop]:= Evaluate(Qech[loop],x[i],V[1][i-omega+1]);
  end for;
end for;

printf "Reduceded System= \n %o \n", Qech;

sol:=[];  // Find values of x_1, \dots, x_{omega-1}
for i:=1 to omega-1 do
  sol[i]:=SquareRoot(MonomialCoefficient(Qech[i],1)/MonomialCoefficient(Qech[i],x[i]^2));
end for;

printf "sol1=%o\n",sol;

// Step 7
for i:=1 to m- omega+1  do
  sol:=sol cat[V[1][i]]; // Append solution of the quadratic system
end for;
 
sol:=sol cat Eltseq(sollin);   // Append solution of the linear equations

printf "sol3=%o\n",sol;

printf "solution in the transformed variables= %o \n", sol;

sol:=VectorSpace( GF,n)!sol;

zet:=sol * Transpose(TMat);  // Apply the Transformation T to find a solution of the original quadratic system

printf "\nSolution of the original system: %o , \n", zet;

for loop := 1 to m do
   ans:= orig[loop] ;
   for i  := 1 to n do
     ans:= Evaluate(ans,x[i],zet[i]);
   end for ;
   printf "ans[%o]=%o\n",loop,ans ;
end for ; 
printf "if all ans=0 then solution satisfies original equations\n";
  
//UnsetOutputFile();