N := 6;  // number of variables
q := 16;   // cardinality of the field

K:=FiniteField(q);
V := VectorSpace(K, N);
R := PolynomialRing(K, N, "grevlex");

printf "N = %o, q = %o\n", N, q;

SetVerbose("Groebner", 2);
SetVerbose("SparseMatrix", 0);


Encrypt := function(PK, x)
    return Vector([Evaluate(p, ElementToSequence(x)) : p in PK]);
end function;
  

function RandomDenseCubicPoly(R)
    cub := &+[Random(K)*R.i*R.j*R.k : i,j,k in [1..Rank(R)]];
    quad := &+[Random(K)*R.i*R.j : i,j in [1..Rank(R)]];
    lin := &+[Random(K)*R.i : i in [1..Rank(R)]];
    cst := Random(K);
    return cub; //+quad+lin+cst;
end function;


function CubicProblemGen()
   // real S and T
  S := Random(GL(N, K));
  U := 1;

  v := Vector(R, [R.i : i in [1..N]]);
  vs := ElementToSequence(v*Matrix(R,S));

  a := [ RandomDenseCubicPoly(R) ];

// make b from a
  b := ElementToSequence(Vector(R, [Evaluate(a[i], vs) : i in [1..U]]));

  return S,a,b;
end function;


// --------------------- GO ! ----------------------------

S,a,b := CubicProblemGen();


x := Random(V);
y := x*S;

assert Encrypt(a, x*S) eq Encrypt(b, x);

db := Encrypt(b, x+y*S^-1) - Encrypt(b, x) - Encrypt(b, y*S^-1);
da := Encrypt(a, x*S+y) - Encrypt(a, x*S) - Encrypt(a, y);

assert da eq db;

Rmat := PolynomialRing(K, 2*N^2, "grevlex");

Rtop := PolynomialRing(Rmat, 2*N, "grevlex");
x := Vector(Rtop, [Rtop.i : i in [1..N]]);
y := Vector(Rtop, [Rtop.i : i in [N+1..2*N]]);

Sx := Matrix(Rmat, N, [Rmat.i : i in [1..N^2]]);
Sx_inv := Matrix(Rmat, N, [Rmat.i : i in [N^2+1..2*N^2]]);
  
inv_eqs := ElementToSequence(Sx*Sx_inv - IdentityMatrix(Rmat, N));
guess := [ Sx[1,i] - S[1,i] : i in [1..N]];

Sx := Matrix(Rtop, Sx);
Sx_inv := Matrix(Rtop, Sx_inv);

if (q lt N) then 
  printf "Adding field equations\n";
  field := [Rmat.i^q - Rmat.i : i in [1..2*N^2]]; 
else
  field := [];
end if;


printf "computing equations\n";


db := Encrypt(b, x+y*Sx_inv) - Encrypt(b, x) - Encrypt(b, y*Sx_inv);
da := Encrypt(a, x*Sx+y) - Encrypt(a, x*Sx) - Encrypt(a, y);

diff_eqs :=  Coefficients((da-db)[1]);

//main_eqs := Coefficients((Encrypt(a, x*Matrix(Rtop,Sx)) - Encrypt(b, x))[1]);


SetVerbose("Groebner", 2);
I := Ideal(diff_eqs cat field cat inv_eqs cat [Rmat.1 - 1]);


V := Variety(I);
#V;