clear;

// *****************************************************************************************************
// *** create a random polynomial quadratic system of m  equations in n variables; 
// *** for solving the system by the method of Hashimoto  one needs  n>=m*(m+3)/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  
// **************************************************************************************************** 
q:=4;
GF<w>:=GaloisField(q);
m:=3;
n:=19;

printf "*********************************************************************** \n";
printf "*** Generate Random System of m equations in n>=m*(m+3)/2 variables *** \n";
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;
b0 := n ge m*(m+3)/2 ;
printf "is n>=m*(m+3)/2: %o\n",b0 ;   

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 "P= %o \n \n", pol;

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

clear ;
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;


printf "*** Parameters ***\n\n";
printf "GF:=GF(%o) \n\n", q;
printf "m:= %o \n\n", m;
printf "n:= %o \n\n\n",n;
Mat,Lin,Con:=PoltoMat(m,n,pol); // build the upper triangular matrices

if n lt (m+3)*m/2 then 
  printf "For m=%o n should be >= %o\n\n",m,m*(m+3)/2;
end if ;


origpol := pol;
origLin := Lin ;
origCon :=Con ;

S:=[];
T:=[];

//printf "original polynomial: %o\n\n",pol;

ell := 1; //----------------------------------------------------------------main loop-----------------------------------
while ell lt m do
  printf "Begin Step %o a : Make the coefficient of x_1x_2 in F1 - Fm to be 0 =================================================================\n\n", ell;
 
  // Ensure that the coefficient of x[ell]^2 <> 0  in F_(m-ell+1)

  S[ell]:=IdentityMatrix(GF,m);

  h:=m-ell+1;
  h0 := h ;
  while Mat[h0][ell][ell] eq 0 do 
    h0 := h0-1;
    if h0 eq 0 then break ;
    end if ;
  end while ;
  if h0 ne 0 then 
    if h ne h0 then 
      temp:= Mat[h0];
      Mat[h0]:=Mat[h];
      Mat[h]:=temp;
      temp := Lin[h0] ;
      Lin[h0] := Lin[h] ;
      Lin[h] := temp ;
      temp := Con[h0] ;
      Con[h0] := Con[h];
      Con[h] := temp ;
      S[ell][h][h]:=0;
      S[ell][h0][h0]:=0;
      S[ell][h0][h]:=1;
      S[ell][h][h0]:=1;
    end if ;

    // Make the coefficient of x1^ell in F_1 ... F_m-ell-1 to 0 
    
    for i:=1 to h-1 do        // Make the coefficient 0
      c:=Mat[i][ell][ell] / Mat[h][ell][ell];
      Mat[i] := Mat[i] - c*Mat[h];
      Lin[i] := Lin[i] - c*Lin[h];
      Con[i] := Con[i] - c*Con[h];
      S[ell][i][h0]:=-c;
    end for;
  else
    printf "For ell=%o the Diagonal term is ZERO !!!!!!!!!!!!!!!!!\n\n",ell;
    printf "matrices:= \n%o\n\n",Mat;
  end if ;


  ell := ell+1 ;


// Make the coefficients Mat[1][ell,1] to Mat[m][ell,ell] to zero

  Tp:=IdentityMatrix(Pol,n);  // Matrix Tp containing unknowns in the ell-th column
  for i:=1 to n do
    Tp[i][ell]:=x[i];
  end for ;


  numb:= 0;
  A1 := [] ;     //ZeroMatrix(GF,numb,n);

  numb1:=0;
  numb2:=0;
  for i:=1 to m+1-ell do       
    for j:=1 to ell-1 do 
      numb1 +:= 1 ;
      for  loop := 1 to n do
        Append(~A1,Mat[i][loop][j]+Mat[i][j][loop]);
      end for ;
    end for ;
  end for ;
  for i:= m+2-ell to m do
    for j:=1 to  m+1-i do 
      numb2 +:= 1 ;
      for loop := 1 to n do 
        Append(~A1,Mat[i][loop][j]+Mat[i][j][loop]);
      end for ;
    end for ;
  end for ;
  numb:=numb1+numb2;
  V:=VectorSpace(GF,numb);
  A1 := Matrix(GF,numb,n,A1)          ;
  
  tn:= NullspaceMatrix(Transpose(A1));
  numb := NumberOfRows(tn) ;  
  // select where t[ell] ne 0
  k := 1 ;
  while tn[k][ell] eq 0 do
    k +:= 1;
    if k gt numb then 
      k-:=1;
      printf "no solution to homogeneous equations %o\n\nNullspace is:%o\n\n",A1,tn;
      break ;
    end if ;
  end while ;
  
  t:=tn[k];

  T[ell]:=IdentityMatrix(GF,n); // Build transformation matrix
  for i:=1 to n do 
    T[ell][i][ell]:=t[i];
  end for;

  for loop:=1 to m do
    Mat[loop]:=Transpose(T[ell]) * Mat[loop] * T[ell]; // Compute new matrix
    // make it upper triangular
    for i :=1 to n do 
       for j:= 1 to i-1 do 
         Mat[loop][j][i] +:= Mat[loop][i][j];
         Mat[loop][i][j] := 0 ;
       end for ; 
    end for ;
    Lin[loop]:=Transpose(T[ell])*Lin[loop];
  end for;

end while;

if Mat[1][m][m] eq 0 then
  printf "Mat[1][%o][%o] is ZERO !!!!!!!!!!!!!\n\n",m,m ;
end if ;

// end of main loop -------------------------------------------------------------------------

pol := MattoPol(Mat,Lin,Con,m,n);
printf "polynomial after step %o :\n%o\n\n",ell,pol;
printf "Transformed matrices for quadratic terms:\n%o\n\n",Mat;

MT:=IdentityMatrix(GF,n);

for i:=2 to m do
  MT:=MT * T[i];
end for;
printf "MT=\n%o\n\n",MT;


MS:=IdentityMatrix(GF,m);

for i:=1 to m-1 do
  MS:=MS * S[i];
end for;
printf "MS=\n%o\n\n",MS;

// Solving the System -----------------------------------------------------

Q:=pol;

Q2:=Q;

Vn:=VectorSpace(GF,m);
we:= Vn!origCon;


// Invert Q
// linear system in x_{m+1}, \dots, x_n
A:=[];
vec:=[];
for loop:=1 to m do
    for i:=1 to m-loop+1 do
       printf  "i:=======%o,  Loop=%o\n",i,loop;
        
        for j:=m+1 to n do
            A:=A cat [MonomialCoefficient(Q[loop],x[i]*x[j])];
            printf " %o x_%o +", MonomialCoefficient(Q[loop],x[i]*x[j]), j;
     
        end for; 
        vec:=vec cat [MonomialCoefficient(Q[loop],x[i])];
        printf " %o =0 \n", MonomialCoefficient(Q[loop],x[i]);
    end for;
end for;

ml:= Integers()!(m*(m+1)/2);
vec:=VectorSpace(GF,ml)!vec;
A0:=Matrix(GF,ml,n-m,A);
b0,t:= IsConsistent(Transpose(A0),vec);
if not assigned t then
  t:= false ;
end if ;

if not b0 then 
    printf "Linear system had no solution matrix A0:\n  %o\nvec=%o\n",Transpose(A0),vec;
else

for b0 := 1 to 1 do
 
  for loop:=1 to m do   // substitute the values of x_m+1, \dots, x_n into Q
    for i:=m+1 to n do
      Q[loop]:=Evaluate(Q[loop],x[i],t[i-m]);
    end for;
  end for;

  printf "\nsystem of form 8.10 Q= %o \n \n", Q; // system of form 8.9

  Q0 :=Q;

  qpol<xx>:=PolynomialRing(GF);
  /* solve  system of form 8.10 from top to bottom */

  sol:=[];
  yps:=[]; // Solve the system 8.9
  for i:=1 to m do
    qp:=MonomialCoefficient(Q[i],x[m-i+1]^2)*xx^2+MonomialCoefficient(Q[i],x[m-i+1])*xx+MonomialCoefficient(Q[i],1);
    solqp := Roots(qp);
    if solqp eq [] then 
      printf "Could not solve system iteratively at step %o\n",i;
      break b0 ;
    end if ;
    sol[i] := solqp[1][1] ;
    // printf "Q[%o]=%o    sol = %o\n\n",i,qp,sol[i];
    for j:=i+1 to m do
      Q[j]:=Evaluate(Q[j],x[m-i+1],sol[i]);
      //printf "Qnew[%o]=%o\n\n",j,Q[j];
    end for;

  end for;
  sol:= Reverse(sol);

  for i:=1 to n-m do // Append the values of x_m+1, \dots, x_n
    sol:=sol cat [t[i]]; 
  end for;

  yps:=VectorSpace(GF,n)!sol;
  zet:= yps*Transpose(MT);  // Invert T

  printf "zet:= %o \n \n", zet;

  printf "Check that we have the correct solution in the new variables \n\n";
  for loop := 1 to m do
    ans:= Q2[loop]  ;
    for i  := 1 to n do
      ans:= Evaluate(ans,x[i],yps[i]);
    end for ;
    printf "ans[%o]=%o\n",loop,ans ;
  end for ;


  printf "Check that we have the correct solution in the original variables \n\n";

  for loop := 1 to m do
    ans:= origpol[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 ;
end for;
end if ;