# required package
using Oscar;

# NxN sudoku
N = 2; F = GF(5); # F = QQ
N = 3; F = GF(11); # F = QQ

Q,T = polynomial_ring(F,"T"); # helper ring
p = prod(T-i for i in (1:N^2))

R,X = polynomial_ring(F,:X => (1:N^2,1:N^2)); # N^2 variables
G = [] # general part

for i in (1:N^2) for j in (1:N^2) # entries in 1:N^2
    G = [ G ; p(X[i,j]) ] end end 

for i in (1:N^2) for j in (1:N^2) for k in (j+1 : N^2) # row conditions
    x1 = X[i,j]; x2 = X[i,k]; G = [ G ; (p(x1)-p(x2))/(x1-x2)] end end end 

for j in (1:N^2) for i in (1:N^2) for k in (i+1 : N^2) # column conditions
    x1 = X[i,j]; x2 = X[k,j]; G = [ G ; (p(x1)-p(x2))/(x1-x2)] end end end

for i in (1:N) for j in (1:N) # box conditions
    for k in (1:N).+(i-1)*N for kk in (1:N).+(i-1)*N 
        for l in (1:N).+(j-1)*N for ll in (1:N).+(j-1)*N  
            if k < kk && l != ll x1 = X[k,l]; x2 = X[kk,ll]; 
               G = [ G ; (p(x1)-p(x2))/(x1-x2)] end end end end end end end

# part depending on M
M = ...
length(findall(x-> x != 0 ,M))

E = G; for i in (1:N^2) for j in (1:N^2) # given entries
           M[i,j] != 0 ?  E = [ E ; X[i,j]-M[i,j] ] : 0 end end 

I = ideal(R, E); # ordering = degrevlex(R)
@time GB = groebner_basis(I, complete_reduction = true);

S = deepcopy(M); rem = []; for i in (1:N^2) for j in (1:N^2)
    nf = normal_form(X[i,j],I); c = constant_coefficient(nf);
    if nf==c S[i,j] = c; else rem = [ rem ; (i,j) ]; 
       print((i,j),": "); nr = N^2*(j-1)+i; d = degree(nf,nr); e = 1;
       if d == 0 println(nf); else while d == e e+=1; 
          nf = normal_form(X[i,j]^e,I); d = degree(nf,nr); end; 
          println(factor(X[i,j]^e-nf)); end end end end; display(S);  

##### examples #####
M = matrix(F,
[ 1 0 0 0;
  3 0 0 0;
  0 0 0 0;
  0 0 0 2; ])
M[3,4] = 1 # = 4

# unique solution (27 entries)
M = matrix(F,
[ 0 6 0 1 0 4 0 5 0;
  0 0 8 3 0 5 6 0 0;
  2 0 0 0 0 0 0 0 1;
  8 0 0 4 0 0 0 0 0;
  0 0 6 0 0 0 3 0 0;
  7 0 0 9 0 1 0 0 4;
  5 0 0 0 0 0 0 0 2;
  0 0 0 2 0 6 9 0 0;
  0 4 0 5 0 8 0 7 0; ])

# leave out an entry, two solutions (26 entries)
M[6,1] = 0 

# leave out an entry (26 entries, takes `forever')
M[8,6] = 0 

# unique solution (26 entries, takes `forever')
M = matrix(F,
[ 0 0 0 0 5 0 0 8 0;
  0 0 0 0 6 2 0 0 5;
  6 0 0 4 0 0 7 0 0;
  0 0 7 0 0 0 9 6 0;
  0 0 5 2 0 6 1 0 0;
  0 3 6 0 0 0 4 0 0;
  0 0 3 0 0 7 0 0 4;
  1 0 0 5 8 0 0 0 0;
  0 6 0 0 1 0 0 0 0; ])
# fix the entries in positions [8,3] and [9,3] 
M[8,3] = 4 # left from 2,4,9 
M[9,3] = 8 # left from 2,8,9