///////// Command:  CII(Input) ///////////////
////// This computes the computable Inoue invariant (CII) of the Sudoku (9 times////// 9) puzzle   

//////////////////// Setting for 9 times 9 Sudoku ///////////////////////////


M2<e1,e2,e3,e4,e5,e6,e7,e8,e9> := PolynomialRing(GF(2),9);
I2:= ideal<M2| e1^2-e1,e2^2-e2,e3^2-e3,e4^2-e4,e5^2-e5,e6^2-e6,e7^2-e7,e8^2-e8,e9^2-e9,e1*e2,e1*e3,e1*e4,e1*e5,e1*e6,e1*e7,e1*e8,e1*e9,
e2*e3,e2*e4,e2*e5,e2*e6,e2*e7,e2*e8,e2*e9,e3*e4,e3*e5,e3*e6,e3*e7,e3*e8,e3*e9,
e4*e5,e4*e6,e4*e7,e4*e8,e4*e9,e5*e6,e5*e7,e5*e8,e5*e9,e6*e7,e6*e8,e6*e9,e7*e8,e7*e9,e8*e9>;
M<e1,e2,e3,e4,e5,e6,e7,e8,e9> := M2/I2;

P<a99,a98,a97,a96,a95,a94,a93,a92,a91,a89,a88,a87,a86,a85,a84,a83,a82,a81,a79,a78,a77,a76,a75,a74,a73,a72,a71,
a69,a68,a67,a66,a65,a64,a63,a62,a61,a59,a58,a57,a56,a55,a54,a53,a52,a51,a49,a48,a47,a46,a45,a44,a43,a42,a41,
a39,a38,a37,a36,a35,a34,a33,a32,a31,a29,a28,a27,a26,a25,a24,a23,a22,a21,a19,a18,a17,a16,a15,a14,a13,a12,a11>:=PolynomialRing(M,81);


/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// Compute the SBG basis of an ideal /////////////////
////////////////////////////////////////////////////////////////////////////

forward Seibun, Boolsudoku, SBGroebner;

SBG2:= function(F)

S := Seibun(F);
G := Boolsudoku(S);
G := SBGroebner(G);

return(Sort(G));

end function;

///////////////////////////////////////////////////////////

Seibun:= function(F)

C := [[],[],[],[],[],[],[],[],[]];
X := [e1,e2,e3,e4,e5,e6,e7,e8,e9];

for t in [1..9] do
	S1 := [m*(P!X[t]) : m in F];

	A := [Coefficients(f) : f in S1| f ne 0];
	B := [Monomials(f) : f in S1 | f ne 0];
	A1:=[];
	    for i in [1..#A] do
		A1 := A1 cat [[Evaluate(m,X[t],1):m in A[i]]];
	    end for;

	for i in [1..#A1] do
		C[t] := C[t] cat [&+[P!(A1[i][j]*B[i][j]):j in [1..#A1[i]]]];
	end for;
end for;

P2<a99,a98,a97,a96,a95,a94,a93,a92,a91,a89,a88,a87,a86,a85,a84,a83,a82,a81,
a79,a78,a77,a76,a75,a74,a73,a72,a71,a69,a68,a67,a66,a65,a64,a63,a62,a61,
a59,a58,a57,a56,a55,a54,a53,a52,a51,a49,a48,a47,a46,a45,a44,a43,a42,a41,
a39,a38,a37,a36,a35,a34,a33,a32,a31,a29,a28,a27,a26,a25,a24,a23,a22,a21,
a19,a18,a17,a16,a15,a14,a13,a12,a11>:=BooleanPolynomialRing(81);

C1:= [P2 ! m : m in C[1]]; 
C2:= [P2 ! m : m in C[2]]; 
C3:= [P2 ! m : m in C[3]]; 
C4:= [P2 ! m : m in C[4]]; 
C5:= [P2 ! m : m in C[5]];
C6:= [P2 ! m : m in C[6]];
C7:= [P2 ! m : m in C[7]];
C8:= [P2 ! m : m in C[8]];
C9:= [P2 ! m : m in C[9]];


C := [C1,C2,C3,C4,C5,C6,C7,C8,C9];

return C;

end function;


////////////////////////////////////////////////////////

Boolsudoku := function(F);

G1 := GroebnerBasis(F[1]);
G2 := GroebnerBasis(F[2]);
G3 := GroebnerBasis(F[3]);
G4 := GroebnerBasis(F[4]);
G5 := GroebnerBasis(F[5]);
G6 := GroebnerBasis(F[6]);
G7 := GroebnerBasis(F[7]);
G8 := GroebnerBasis(F[8]);
G9 := GroebnerBasis(F[9]);



G11 := [e1*(P!m) : m in G1];
G12 := [e2*(P!m) : m in G2];
G13 := [e3*(P!m) : m in G3];
G14 := [e4*(P!m) : m in G4];
G15 := [e5*(P!m) : m in G5];
G16 := [e6*(P!m) : m in G6];
G17 := [e7*(P!m) : m in G7];
G18 := [e8*(P!m) : m in G8];
G19 := [e9*(P!m) : m in G9];


G20 := G11 cat G12 cat G13 cat G14 cat G15 cat G16 cat G17 cat G18 cat G19;

return G20;
end function;

//////////////////////////////////////////////

SBGroebner := function(G) ///// stratification

X := {@ LeadingMonomial(m) : m in G @};
n := #X;
Y := [i*(P ! 0) : i in [1..n]];

for i in [1..n] do
    Y[i] := &+[x : x in G | LeadingMonomial(x) eq X[i]];
end for;

return Y;
end function;

///////////////////////////////////////////////

Const := function(F);
Zero := [i*0:i in [1..81]];
U := [m: m in F | Exponents(LeadingTerm(m)) eq Zero];
return #U; 
end function;


//////////////////////////////////////////////////////
////////////////////////////////////////////////////
/////////極小多項式関係

/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
////////////////////////////////////////////////////
////// 多項式fを入力して，係数列のwedge積を出力する．
/////////////////////////////////////////////////

wed := function(a,b);
d:= a+b+a*b;
return d;
end function;

wedge:= function(f);
F:= Coefficients(f);
s:= #(F);
A:= F[1];
if s eq 1 then return A;
else 
for i in [2..s] do
A:= wed(A,F[i]);
end for;
end if;

return A;

end function;

///////////////////////////多項式fから定数項を取り出すコマンド．
/////////////////定数項が0なら0を返す．

teisuu := function(f);
t:= Rank(Parent(f));
Z:= [i*0: i in [1..t]];
p:= [r: r in Terms(f) | Exponents(r) eq Z];

if p eq [] then return Parent(f)!0;
else return p[1];
end if;

end function;



////////////////////////////////////////////////////////////////
//// 多項式の配列 F から，１次式のみを抜き出して，それらのminipoly の配列を出力
/////////////////////////////////////////////////////////////

minipoly := function(F);
F:= [f: f in F | LeadingTotalDegree(f) eq 1];
s:= #(F);
E:= e1+e2+e3+e4+e5+e6+e7+e8+e9;

for i in [1..s] do 
g:= F[i] - LeadingTerm(F[i]) - teisuu(F[i]);

if g eq 0 then F[i] := F[i];
else 
h:= LeadingCoefficient(F[i])*(E+wedge(g))*LeadingMonomial(F[i])+ teisuu(F[i])*(E+wedge(g)); 
F[i]:= h;
end if;

end for;

return F;

end function;


/////////////////////////////////////////////////////////////
/////////// 

///// 極小多項式の配列から2種類のcontradiction 多項式全部を取り出す命令

contra:= function(F);

F11 := [f: f in F| teisuu(f) eq 0 and Length(LeadingCoefficient(f)) eq 9];
F21 := [f: f in F | teisuu(f) ne 0 and Length(teisuu(f)) ge 2];

G:= F11 cat F21;

return G;

end function;

///////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////　極小多項式の配列からstrategy多項式の配列を取り出す命令

refuta := function(F);
F1 := [f: f in F | teisuu(f) eq 0 and Length(LeadingCoefficient(f)) le 7];

return F1;

end function;

///////////////////////////////////////////////////////////
//////strategy多項式の配列から，その配列に含まれる先頭係数の個数の配列の和を出力

sentou:= function(F);

if F eq []
then return 0;
end if;

F1 := [Length(LeadingCoefficient(f)): f in F];

return &+(F1);

end function;

/////////////////////////////////////////////////////////
////極小多項式の配列から、まだ値の定まっていない変数の配列を出力

unvariable := function(F);

F1 := refuta(F);
F2:= [LeadingMonomial(f): f in F1];

return F2;

end function;

////////////////////////////////////////////////////////
/////// 極小多項式の配列から，semi-sol 多項式の配列を取り出す命令

semipoly := function(F);

E:= e1+e2+e3+e4+e5+e6+e7+e8+e9;

F1:= [f: f in F | teisuu(f) ne 0 and Length(teisuu(f)) eq 1 and LeadingCoefficient(f) ne E];
F11 := [f: f in F1 | teisuu(f)*LeadingCoefficient(f) ne 0];

F2 := [f: f in F | teisuu(f) eq 0];
F21 := [f: f in F | Length(LeadingCoefficient(f)) eq 8];

G := F11 cat F21;
G := SetToSequence(SequenceToSet(G));

return Sort(G);

end function;

/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
////////////極小多項式の配列から，付け加えるべきsol. poly の配列を出力

solpoly := function(F);


E:= e1+e2+e3+e4+e5+e6+e7+e8+e9;

F1:= [f: f in F | teisuu(f) ne 0 and Length(teisuu(f)) eq 1 and LeadingCoefficient(f) ne E];
F11 := [f: f in F1 | teisuu(f)*LeadingCoefficient(f) ne 0];
F12 := [E*LeadingMonomial(f) + teisuu(f): f in F11];

F2 := [f: f in F | teisuu(f) eq 0];
F21 := [f: f in F2 | Length(LeadingCoefficient(f)) eq 8];
F22:= [E*LeadingMonomial(f) +(E-LeadingCoefficient(f)) : f in F21];

G:= F12 cat F22;
G := SetToSequence(SequenceToSet(G));

return G;

end function;

///////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
////極小多項式の配列から，sol. poly のみを取り出して配列を作る．

solutionpoly := function(F);

E:= e1+e2+e3+e4+e5+e6+e7+e8+e9;
F1:= [f: f in F | teisuu(f) ne 0];
F2 := [f: f in F1| LeadingCoefficient(f) eq E and Length(teisuu(f)) eq 1];

return F2;

end function;



//////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////
///// 極小多項式の配列から、semipoly とsolpoly の合計数を出力
AlmostSolution4 := function(F); //// output #(SP(F) \cup Sol(ASP(F))) 

aa:= minipoly(F);
bb:= #(semipoly(aa) cat solutionpoly(aa));

return bb;
                              
end function;

////////////////////////////////////////////////////////////////////////////////

 
/////////////////////////////////////////////////////////
//////////////// SBG基底を入力して，sol. poly達を付け加えて，次のSBG基底を出力
//////// １段階だけ////見事に付け加えていくsemipoly の数は旧SBGTransformと同じ！///// minipoly を使っているので、G に定数がはいっていないときのみ使用可

SBGTrans2 := function(G);

H := solpoly(minipoly(G));
H1 := SBG2(G cat H);

return H1;

end function;

///////////////////////////////////////////
/////////basic closure を出力するプログラム

SBGTrans := function(F);

G := SBG2(F);
if Const(G) ne 0 then return G;
else mini := minipoly(G);
    if contra(mini) ne [] then return G;
    else 

while semipoly(mini) ne [] do 

G := SBGTrans2(G);
if Const(G) ne 0 then return G;
else mini:= minipoly(G);
     if contra(mini) ne [] then return G;
     end if;
end if;
     
end while;

return G;

end if;
end if;

end function;

//////////////////////////////////////////////////////



/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
////////////////////////////////////////////////////
forward Solpoly, variable,Const, nakano8;
nakano8 := function(F: d := 0, W := []) 
/// The main body of this program /////////////

L := [* *];

X := [e1,e2,e3,e4,e5,e6,e7,e8,e9];
E := e1+e2+e3+e4+e5+e6+e7+e8+e9;
Zero := [i*0:i in [1..81]];

G:= SBGTrans(F);

S:= Solpoly(G); 
"#(solution polynomials) = ", #S;

d:= d+1;
"level=", d;


if Const(G) ne 0 then L := L cat [*[* 1,d,W *]*];
print "NSleaf_1","//////////";
else if contra(minipoly(G)) ne [] then L := L cat [*[* 1,d,W *]*];
        print "NSleaf_2", "//////////";
     else if  #solutionpoly(minipoly(G)) eq 81 then L := L cat [*[* 1,d,W *]*];
           print "SolutionLeaf","//////////";
          else 

mini := minipoly(G);

V := unvariable(mini); ///the sequence of undetermined variables
B := [];
H := [];
H3 := [];

/////////////////////////////
for v in V do 
Z:= [f: f in mini | LeadingMonomial(f) eq v];
z:= Z[1];
X1:= [x: x in X | LeadingCoefficient(z)*x eq 0];

K :=[G cat [E*v+x]: x in X1];
K1:= [SBG2(K[i]): i in [1..#X1]];
/// K2:= [m : m in K1 | Const(m) eq 0 or contra(minipoly(m)) ne []];
B := B cat [#X1];

///H1:= [i: i in [1..#X1] | K1[i] in K2];
///H2:= [X1[i]: i in H1];
H := H cat [X1]; 
H3:= H3 cat [K1];
end for;
/// V;B;H; 
////////////////////////////////////

b := Min(B);///the minimal number of branches

if b eq 0 then L := L cat [*[* 1,d,W *]*];// one more leaf
print "NSleaf_3","//////////";/// No solution leaf
end if;

if b ne 0 then 

V1 := [j: j in [1..#V] | B[j] eq b]; ///分岐数が最小（b個）の番号集合
V2 := [V[j]: j in V1]; ///分岐数が最小の変数の集合
H4 := [H[j]: j in V1]; /// 各変数について取りうる値の配列を並べた集合
H5 := [H3[j]: j in V1]; /// 各変数について取りうる値を加えたSBG基底の配列の集合
H6:= [[AlmostSolution4(H5[i][j]): j in [1..b]] : i in [1..#V1]];
  
H7:= [&+(H6[i]): i in [1..#V1]];

b1 := Max(H7);
H9 := [i: i in [1..#V1] | &+(H6[i]) eq b1]; // 和が最大の番号の配列

V2 := [V2[j]: j in H9]; 
H4 := [H4[j]: j in H9]; 
H5 := [H5[j]: j in H9]; 


H8:= [[minipoly(H5[i][j]): j in [1..b]] : i in [1..#H9]];
H81 := [[refuta(H8[i][j]): j in [1..b]] : i in [1..#H9]];
H82 := [[sentou(H81[i][j]): j in [1..b]] : i in [1..#H9]];
H83:= [&+(H82[j]): j in [1..#H9]];

b2:= Max(H83);
H11:= [j: j in [1..#H9] | H83[j] eq b2];

H12 := [V2[j]: j in H11]; /// 最適の変数集合
H14:= [H4[j]: j in H11]; /// [e2,e4,e5] などの配列
H15:= [H5[j]: j in H11]; 

H13 := [];
for i in [1..#H11] do 
H13 := H13 cat [H15[i]];
end for;

L := L cat [*[* 0,d,W *]*];

print "[#(branches),#(variables)]=",[b,#H11], "/////";

for j in [1..#H12] do 
for k in [1..b] do
L1:= nakano8(H13[j][k]: d:= d, W := [[ H12[j], H14[j][k] ]] cat W); 
L := L cat L1;
end for;
end for;

end if;

end if; //
end if;
end if;

return L;

end function;
//////////////////////////////////////////////////////
///////////////////////////////////////////////////
////////////////////////////////////////////////////
/////////// 井上樹形図

ITR:= function(L)

H1:= [x[2]: x in L];
d:= Max(H1);

if d eq 1 then return [[2,1,1]];
end if;

H1:= [];
for k in [1..d] do
H2:= [x[3]: x in L | x[2] eq k];
H1:= H1 cat [H2];
end for;


S := [];
K:= H1[2];
K2:= {x[1][1]: x in K};
k2:= #K2;
for z in K2 do 
K3:= [[]] cat [y: y in K| y[1][1] eq z];
S:= S cat [K3];
end for;
if d eq 2 then S:= S;
end if;

for i in [3..d] do
K:= H1[i];


for x in K do 

S1:= [y: y in S | Remove(x,1) in y];

//"S1=",S1;
S2:= [y: y in K | Remove(y,1) eq Remove(x,1)];

//"S2=",S2;
S3:= {z[1][1]: z in S2}; 

for w in S3 do
S4:= [x: x in S2 | x[1][1] eq w];


for z in S1 do 
S:= S cat [z cat S4];
S:= [x: x in S | x notin S1];
end for;
end for;


K:= [x: x in K | x notin S2];
if K eq [] then break;
end if;
end for;

end for;

e:= #S;

if e eq 1 then 
H10 := [[]];
end if;

H10:= [[]];
for i in [1..e-1] do
H10 := H10 cat [[]];
end for;

for i in [1..e] do 
H10[i] := [x: x in L | x[3] in S[i]];
end for;

K5:= [];
for i in [1..e] do
K6:= [x: x in H10[i] | x[1] eq 1]; 
l:= Max([x[2]: x in H10[i]]);
K5:= K5 cat [[#H10[i]+1,#K6,l]];
end for;

//H10;
return K5;

end function;

//////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////




////////////////////////////////////////////////

Solpoly := function(F); ///output the set of solution polynomials 
                        ///contained in the SBG bases F

X := [e1,e2,e3,e4,e5,e6,e7,e8,e9];
E := e1+e2+e3+e4+e5+e6+e7+e8+e9;
Y := [E*a11,E*a12,E*a13,E*a14,E*a15,E*a16,E*a17,E*a18,E*a19,E*a21,E*a22,E*a23,E*a24,E*a25,E*a26,E*a27,E*a28,E*a29,
      E*a31,E*a32,E*a33,E*a34,E*a35,E*a36,E*a37,E*a38,E*a39,E*a41,E*a42,E*a43,E*a44,E*a45,E*a46,E*a47,E*a48,E*a49,
      E*a51,E*a52,E*a53,E*a54,E*a55,E*a56,E*a57,E*a58,E*a59,E*a61,E*a62,E*a63,E*a64,E*a65,E*a66,E*a67,E*a68,E*a69,
      E*a71,E*a72,E*a73,E*a74,E*a75,E*a76,E*a77,E*a78,E*a79,E*a81,E*a82,E*a83,E*a84,E*a85,E*a86,E*a87,E*a88,E*a89,
      E*a91,E*a92,E*a93,E*a94,E*a95,E*a96,E*a97,E*a98,E*a99];

E := e1+e2+e3+e4+e5+e6+e7+e8+e9;
S := [m : m in F | LeadingTerm(m) in Y and #Terms(m) eq 2 and TrailingTerm(m) in X];
     // the set of solution polynomials

return S;

end function; 

/////////////////////////////////

variable := function(F);
Y := [a11,a12,a13,a14,a15,a16,a17,a18,a19,a21,a22,a23,a24,a25,a26,a27,a28,a29,
      a31,a32,a33,a34,a35,a36,a37,a38,a39,a41,a42,a43,a44,a45,a46,a47,a48,a49,
      a51,a52,a53,a54,a55,a56,a57,a58,a59,a61,a62,a63,a64,a65,a66,a67,a68,a69,
      a71,a72,a73,a74,a75,a76,a77,a78,a79,a81,a82,a83,a84,a85,a86,a87,a88,a89,
      a91,a92,a93,a94,a95,a96,a97,a98,a99];

S := Solpoly(F);
K := [LeadingMonomial(m) : m in S]; 
K1 := [m : m in Y | m notin K];
return K1;

end function;


///////////////////////////////////////////
Const := function(F);
Zero := [i*0:i in [1..81]];
U := [m: m in F | Exponents(LeadingTerm(m)) eq Zero];
return #U; 
end function;


////////////////////////////////
///////////////////////////////
////// CII を計算するプログラム
///////////////////////////////

CII := function(F)

F:= SBG2(F);
L:= nakano8(F);
N := ITR(L);
return Min(N), #N;

end function;

//////////////////////////////////////
//////////////////////////////////////