rel:= function(a,b,c,d)
  A<X>:=PolynomialAlgebra(Parent(a));
  return Coefficients((X-a)*(X-b)*(X-c)*(X-d)-(X-1)*(X-2)*(X-3)*(X-4));
end function;

sudokujr := function (cst)
 local A,H,V,C,L,Q,map;
 A:=PolynomialAlgebra(GF(5),16);
 H:=&cat [rel(A.(4*i+1),A.(4*i+2),A.(4*i+3),A.(4*i+4)): i in [0..3]];
 V:=&cat [rel(A.(i+1),A.(i+5),A.(i+9),A.(i+13)): i in [0..3]];
 C:=&cat [rel(A.(2*i+8*j+1),A.(2*i+8*j+2),A.(2*i+8*j+5),A.(2*i+8*j+6)):
           i in [0..1], j in [0..1]];
 E:=[A.i-cst[i]: i in [1..16] | cst[i] ne 0];
 G:=GroebnerBasis(H cat V cat C cat E);
 print("");
 return Matrix(A,4,4,[NormalForm(A.i,G): i in [1..16]]);
end function;

sudokujr([0,0,0,4,\
          0,2,0,0,\
          0,0,3,0,\
          1,0,0,0]);

sudokujr([1,0,0,3,\
          4,0,0,0,\
          0,0,0,0,\
          0,0,0,2]);

sudokujr([0,0,0,0,\
          0,1,2,0,\
          0,2,1,0,\
          0,0,0,0]);

sudokujr([1,0,0,0,\
          0,0,0,4,\
          0,1,0,0,\
          0,0,0,3]);

rel2:= function(L)
  A<X>:=PolynomialAlgebra(Parent(L[1]));
  return Coefficients(&*[X-a:a in L]-&*[X-i:i in [1..#L]]);
end function;

sudoku := function (cst)
 local m,m2,A,S,D,d,H,V,C,L,Q,map;
 A:=PolynomialAlgebra(GF(11),81);
 c:=function(i)
   if cst[i] eq 0 then
     return A.i;
   else 
     return cst[i];
   end if;
 end function;
 H:=&cat [rel2([c(9*i+j):j in [1..9]]): i in [0..8]];
 V:=&cat [rel2([c(i+j*9):j in [0..8]]): i in [1..9]];
 C:=&cat [rel2([c(3*i+27*j+k+9*l):k in [1..3], l in [0..2]]): i in [0..2], j in [0..2]];
 E:=[A.i-cst[i]: i in [1..81] | cst[i] ne 0];
 G:=GroebnerBasis(H cat V cat C) cat E;
 return Matrix(A,9,9,[NormalForm(A.i,G): i in [1..81]]);
end function;
SetVerbose("Groebner",1);

//sudoku([2,5,8,0,3,0,9,0,1,\
        0,1,0,8,0,4,0,0,0,\
        4,0,7,0,0,0,2,0,8,\
        0,0,5,2,0,0,4,8,6,\
        7,0,0,0,9,8,1,0,0,\
        8,4,0,0,0,3,0,0,0,\
        1,0,0,3,6,0,0,7,2,\
        0,7,0,0,0,0,0,0,3,\
        9,0,3,5,8,7,6,0,4]);

//sudoku([0,0,3,0,0,4,6,0,0,\
        2,0,6,0,9,0,0,0,0,\
        7,0,0,0,8,3,0,5,1,\
        0,4,0,8,0,0,7,0,9,\
        5,7,0,3,6,0,1,4,0,\
        9,0,2,0,0,1,5,0,0,\
        0,0,5,9,0,6,0,0,0,\
        0,0,0,4,2,0,0,0,5,\
        0,0,9,5,0,7,4,3,0]);
