//Written by Enrico Thomae. Please note that this is a preliminary and naive implementation. Questions or mistakes can be send to enrico.thomae@rub.de.

//------------------
//----Parameters----
//------------------

m:=9;           // Number of equations
n:=3*m;         // Number of variables
exp:=8;
k:=2^(exp);     // GF(k)

"-----------------";
"m=",m,", w=",w,", k=",k;
"-----------------";

//--------------
//----Macros----
//--------------

PolyRing<[t]> := PolynomialRing(GF(k),n);
PolyBigRing<[s]> := PolynomialRing(GF(k),n^2);
tau := Binomial(n+1,2);  // Number of quadratic monomials

//----Transformation of Variables (Sigma)----

e:=[];
for i in [1..n] do
  for j in [i..n] do
    Append(~e, [i,j]);
  end for;
end for;

function Trafo(T)
  K :=ZeroMatrix(PolyBigRing,tau,tau);
  for p in [1..tau] do
    if e[p][1] eq e[p][2] then
      for i:=1 to tau do
        K[p,i]:=T[e[i,1],e[p][1]]*T[e[i,2],e[p][2]];
      end for;
    else
      for i:=1 to tau do
        K[p,i]:=T[e[i,1],e[p][1]]*T[e[i,2],e[p][2]]+T[e[i,1],e[p][2]]
*T[e[i,2],e[p][1]];
      end for;
    end if;
  end for;
  K := Transpose(K);
  return K;
end function;

//----label(x,y) returns the corresponding row of Sigma----

function label(x,y)
  sum:=0;
  if (x eq 1) then return y; end if;
  if (x gt 1) then 
    for temp:=0 to x-2 do sum+:=(n-temp); end for;
    return sum+y-x+1;
  end if;
end function;

//----pos(x,y) enumerates the entries of the matrix S----

function pos(x,y)
  return n*(y-1)+x;
end function;

//----Test solution----

function testsolution(MQtest,solvec)
  test := true;
  for temp := 1 to #MQtest do
    if Evaluate(MQtest[temp],solvec) ne 0 then test := false; end if;
  end for; 
  return test;
end function;

//----Test transformation----

function testtrafo(MQ,S)
  MQtrafo := [];
  vartraf:=[];

  for temp1 := 1 to n do
    curPoly := Zero(PolyRing);
    for temp2 := 1 to n do
       curPoly+:=GF(k)!S[temp1,temp2]*t[temp2];
    end for;
    Append(~vartraf,curPoly);
  end for;

  for temp1 := 1 to #MQ do
    Append(~MQtrafo,Evaluate(MQ[temp1],vartraf));
  end for;
 
  return MQtrafo;
end function;

//----Find solution of Groebner Basis----

function givesol(B,sol)
  if B[1] eq 1 then err := true;
    else err := false;
  end if;
  for temp := #B to 1 by -1 do
    t1,t2,t3 := IsUnivariate(B[temp]);
    if t1 then 
      s1,s2 := HasRoot(UnivariatePolynomial(B[temp]));
      if s1 then
        sol[t3] := s2;
      else err := true; break temp;
      end if;
    else err := true; break temp;
    end if;
    for temp2 := 1 to #B do
      B[temp2] := Evaluate(B[temp2],sol);
    end for;
  end for;
  return sol,err;
end function;

//--------------------------------------------
//----Initialisation of random MQ instance----
//--------------------------------------------

MQ := [];

for temp := 1 to m do
  curPoly := Zero(PolyRing);
  for i := 1 to n do
    for j := i to n do
      curPoly +:=Random(GF(k))*t[i]*t[j];
    end for;
  end for;
  for i := 1 to n do
    curPoly +:=Random(GF(k))*t[i];
  end for;
  curPoly +:=Random(GF(k));
  Append(~MQ, curPoly);
end for;

MQCoeff := [];

for temp := 1 to m do
  templist := [];
  for i := 1 to n do 
    for j := i to n do
      Append(~templist, MonomialCoefficient(MQ[temp],t[i]*t[j]));
    end for;
  end for;
  Append(~MQCoeff, templist);
end for;

//------------------------------------------------------
//----Standard approach: Guessing and Groebner Basis----
//------------------------------------------------------

" ";
"----Standard approach----";

guess := [];
for temp := 1 to m do Append(~guess,t[temp]); end for;
for temp := m+1 to n do Append(~guess,Random(GF(k))); end for;

MQguess := [];
for temp := 1 to m do Append(~MQguess,Evaluate(MQ[temp],guess)); end for;

reduceSeq := [PolyRing ! t[i]^k - t[i] : i in [1..m]]; // use only if field size is smaller than degree of regularity 
GBideal := ideal<PolyRing |MQguess>;

curTime := Cputime();
B := GroebnerBasis(GBideal);
curTime := Cputime(curTime);

guess,err := givesol(B,guess);

if err then 
  "No solution found. Try again."; 
else "Time to solve with standard approach: ", curTime;
  "Solution: ", guess; // print solution
  "Test of solution: ", testsolution(MQ,guess);
end if;

//-----------------------------------
//----Approach of Thomae and Wolf----
//-----------------------------------

" ";
"----Approach of Thomae and Wolf----";

//----Determine linear transformation S----

S := ScalarMatrix(n,PolyBigRing!1);
for j := 1 to m do
  for i := 1 to n do
    S[i,j] := s[pos(i,j)];
  end for;
end for;

for i := 1 to n do 
  S[i,1] := Random(GF(k));
end for;

w := Truncate(n/m)-1;
err := true;
timeTW := 0;

for curcol := 2 to m do

  if w*(curcol-1) lt n then 
    for temp := 1 to n-w*(curcol-1) do
      S[temp,curcol] := Random(GF(k));
    end for;
  end if;

  Gls:=[];
  for i := 1 to curcol-1 do
    for j:=1 to w do
      tempGl := Zero(PolyBigRing);
      tempMa := Trafo(S);
      for temp:=1 to tau do 
        tempGl+:=MQCoeff[j,temp]*tempMa[temp, label(i,curcol)]; 
      end for;
      Append(~Gls, tempGl);
    end for;
  end for;

  GBideal := ideal<PolyBigRing | Gls>;
  curTime := Cputime();
  B := GroebnerBasis(GBideal);
  curTime := Cputime(curTime);
  timeTW +:= curTime;

  err := false;
  for temp := 1 to #B do
    if (#Terms(B[temp]) gt 2) or (B[1] eq 1) then err := true; end if; 
  end for;

  if err then
    "Linear system of equations contains no or more than one solution. Try again."; break curcol;
  else
    for temp := 1 to #B do
      S[n-w*(curcol-1)+temp,curcol] := -MonomialCoefficient(B[temp],1);
    end for;
  end if;

end for;

//----Determine vinegar variables----

if not(err) then

  MQtrafo := testtrafo(MQ,S);
  Gls := [];

  for temp1 := 1 to w do
    for temp2 := 1 to m do
      tempGl := Zero(PolyRing);
      for temp3 := m+1 to n do 
        tempGl +:= MonomialCoefficient(MQtrafo[temp1],t[temp2]*t[temp3])*t[temp3]; 
      end for;
      tempGl +:= MonomialCoefficient(MQtrafo[temp1],t[temp2]);
      Append(~Gls, tempGl);
    end for;
  end for;

  GBideal := ideal<PolyRing | Gls>;
  curTime := Cputime();
  B := GroebnerBasis(GBideal);
  curTime := Cputime(curTime);
  timeTW +:= curTime;

  err2 := false;
  for temp := 1 to #B do
    if (#Terms(B[temp]) gt 2) or (B[1] eq 1) then err2 := true; end if; 
  end for;

  if err2 then 
    "No or more than one solution found. Try again."; 

    else 

    var:=[];
    for temp := 1 to m do
      Append(~var,t[temp]);
    end for;
    for temp := m+1 to n do
      Append(~var,-MonomialCoefficient(B[temp-m],1));
    end for;
 
    for temp := 1 to m do
      MQtrafo[temp] := Evaluate(MQtrafo[temp],var);
    end for; 

    //----Obtain w linear equations----

    for temp := 1 to w do

      tempPoly := MQtrafo[temp]^(2^(exp-1));

      tempGl := Zero(PolyRing);
      for temp2 := 1 to m do 
        tempGl +:= MonomialCoefficient(tempPoly,t[temp2]^k)*t[temp2]; 
      end for;
      tempGl +:= MonomialCoefficient(tempPoly,1);

      MQtrafo[temp] := tempGl;
    end for;

    GBideal := ideal<PolyRing | MQtrafo>;
    curTime := Cputime();
    B := GroebnerBasis(GBideal);
    curTime := Cputime(curTime);
    timeTW +:= curTime;

    var,err := givesol(B,var);
    
    if err then 
      "No solution found. Try again."; 
    else "Time to solve with TW approach: ", timeTW;
      //----Transform solution back with S
      solution := [];
      for temp3 := 1 to n do
        tempsol := GF(k)!0;
        for temp4 := 1 to n do
          tempsol +:= GF(k)!S[temp3,temp4]*GF(k)!var[temp4];
        end for;
        Append(~solution,tempsol);
      end for;
      "Solution: ", solution; // print solution
      "Test of solution: ", testsolution(MQ,solution);
    end if;   
  end if;
end if;























