/* searching for winning-strategy for Euler Getter */
/* since 2012/02/17 */

/* how to use

load:
batchload("c:/ ..(pass).. /EG.txt");

start:
AutoTurnStart(player);
Here player=+1 if you are red, or 1st player
            -1 if you are blue, or 2nd

continue:
AutoTurn(x,y);
is to take place at (x,y), here (x,y) is like
(x,y)-entry of a matrix.
Upper-right point is also referred as (0,0).

parameters:
To change size of board, you need to change the value of N
by rewriting the line below.
After that, you must reload.
(Technically, you can also change parameters R1,R2,R3
related to Monte-Carlo method.)

*/

/* size of diagram */
/* the number of cells M is N*N+1 */
N:4;
M:N^2+1;

/* data of tableau */
/*
 a b c ...
 d e f ...
 .........
 g h i ...
 z <- A-point

 we represent this tableau by the array
 [a,b,c,...,z]
 entry=+1 if red (or the first player)
       -1 if blue (or the second)
        0 if void
*/

/*******************************/
/* computation of euler number */
/*******************************/

/* we compute euler number by adding cell one-by-one */

/* patching xy-plane */
/* be careful around A&B-points */
NormalizeXY(xy):=block(
 [x,y],
 x:xy[1],y:xy[2],
 if xy=[0,N] or xy=[N,0] then return([N,0]),
 if xy=[N,N] then return([0,0]),
 /* assumption: x or y =1..N-1 */
 if x=-1 then(x:mod(x,N),y:mod(-y-1,N)),
 if x= N then(x:mod(x,N),y:mod(-y,N)),
 if y=-1 then(y:mod(y,N),x:mod(-x-1,N)),
 if y= N then(y:mod(y,N),x:mod(-x,N)),
 return([x,y])
);

/* list of adjacent cells in clockwise order */
/* N>2 */
AdjacentCells(i):=block(
 [x,y,ret,j],
 y:mod(i-1,N),
 x:(i-1-y)/N, /* 0<=x,y<N except for A-point */
 /* A-point */
 if i=N^2+1 then(return([N,N*(N-1)+1])),
 /* B-point */
 if i=1 then(return([2,N+2,N+1,N^2])),
 /* other points */
 ret:[ [x,y+1],[x+1,y+1],[x+1,y],[x,y-1],[x-1,y-1],[x-1,y] ],
 for j:1 thru 6 do(
  ret[j]:NormalizeXY(ret[j]),
  ret[j]:ret[j][1]*N+ret[j][2]+1
 ),
 return(ret)
);

MakeAdjList():=block(
 [AdjList],
 AdjList:[],
 for i:1 thru M do(
  AdjList:append(AdjList,[AdjacentCells(i)])
 ),
 return(AdjList)
);
/* to make AdjList */
AdjList:MakeAdjList();

EulerNumber(t,player):=block(
 [e,t0,i,j,adj],
 e:0, /* euler number */
 t0:makelist(0,i,1,M),
 if N=1 then(
  if t[2]=player then e:e+1,
  return(e)
 ),
 for i:1 thru M do(
  if t[i]=player then(
   adj:AdjList[i],
   for j:1 thru length(adj)-1 do(
    if t0[adj[j]]#player and t0[adj[j+1]]=player then(e:e-1)
   ),
   if t0[adj[length(adj)]]#player and t0[adj[1]]=player then(e:e-1),
   e:e+1, /* cf. Mayer-Vietoris sequence */
   t0[i]:player
  )
 ),
 return(e)
);

Winner(t):=block(
 if EulerNumber(t,+1)>0 then return(+1)
 else(return(-1))
);

/**********************************************/
/* to determine victorious subset inductively */
/**********************************************/

/* to make the set S[n] of states at step n */

SetIncrement(set,max):=block(
 [i,j,r,flag],
 r:length(set),
 flag:false,
 if set[r]<max
 then(
   set[r]:set[r]+1
 )else(
  for i:r-1 thru 1 step -1 do(
   if set[i]<set[i+1]-1 then(
    set[i]:set[i]+1,
    for j:i+1 thru r do(
     set[j]:set[j-1]+1
    ),
    i:-1,
    flag:true
   )
  ),
  if flag=false then(set:makelist(i,i,1,r))
 ),
 return(set)
);

MakeS(n):=block(
 [S,k,set1,set2,t,i,c],
 S:[],
 k:ceiling(n/2),
 c:0,
 set1:makelist(i,i,1,n), /* n-subset of {1..M} */
 do(
  set2:makelist(i,i,1,k), /* k-subset of {1..n} */
  do(
   t:makelist(0,i,1,M),
   for i:1 thru n do(t[set1[i]]:-1),
   for i:1 thru k do(t[set1[set2[i]]]:+1),
   S:append(S,[t]),
   c:c+1,if mod(c,10000)=0 then print("---",c),
   if set2[1]>n-k then return(),
   SetIncrement(set2,n)
  ),
  if set1[1]>M-n then return(),
  SetIncrement(set1,M)
 ),
 return(S)
);

MakeFiltS(fil,n):=block(
 [M0,n0,pos,S,k,set1,set2,t,i,c],
 pos:[],
 for i:1 thru M do(
  if fil[i]=0 then pos:append(pos,[i])
 ),
 M0:length(pos),
 n0:n-(M-M0),
 S:[],
 if mod(M-M0,2)=0 then k:ceiling(n0/2) else k:floor(n0/2),
 c:0,
 set1:makelist(i,i,1,n0), /* n0-subset of {1..M0} */
 do(
  set2:makelist(i,i,1,k), /* k-subset of {1..n} */
  do(
   t:copylist(fil),
   for i:1 thru n0 do(t[pos[set1[i]]]:-1),
   for i:1 thru k do(t[pos[set1[set2[i]]]]:+1),
   S:append(S,[t]),
   c:c+1,if mod(c,10000)=0 then print("---",c),
   if length(set2)=0 or set2[1]>n0-k then return(),
   SetIncrement(set2,n0)
  ),
  if set1[1]>M0-n0 then return(),
  SetIncrement(set1,M0)
 ),
 return(S)
);

/* filtering of result */

Filt(t,U):=block(
 [T,num,i,j,flag],
 T:[],
 num:0,
 for i:1 thru length(U) do(
  flag:true,
  for j:1 thru M do(
   if t[j]#0 and t[j]#U[i][j] then(flag:false,j:M+1)
  ),
  if flag=true then(T:append(T,[copylist(U[i])]))
 ),
 return(T)
);

FiltPrint(t):=block(
 [ct],
 for ct:1 thru M do(
  print("#S_A(",ct,")=",length(Filt(t,Sa[ct])),
        ",#S_B(",ct,")=",length(Filt(t,Sb[ct])) )
 )
);

/* to make victorious subset for the end of game (n=M) */
/* * * * *
print("#S(",M,")=",binomial(M,ceiling(M/2)));
/*S[M]:MakeS(M);*/

fil:[-1,0,0,1, 0,0,0,0, 0,0,0,0, -1,0,0,0, 1];
fil:[1,1,0,0, 0,-1,0,0, 0,0,0,0, 0,0,0,1, -1];
fil:[- 1, 1, 1, 1, 1, 1, - 1, 0, - 1, - 1, - 1, 1, - 1, - 1, 0,
     - 1, - 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, - 1, - 1, 1, - 1,
     0, - 1, 0, - 1, 1];
/*fil:makelist(0,i,1,M);fil[1]:+1;*/
print("filtered by",fil);
PrintTab(fil);
S[M]:MakeFiltS(fil,M);

Sa[M]:[]; /* for the first player */
Sb[M]:[]; /* for the second */
for i:1 thru length(S[M]) do(
 e:EulerNumber(S[M][i],+1),
 if e>0 /*e<=0*/ then(
  Sa[M]:append(Sa[M],[S[M][i]])
 )else(
  Sb[M]:append(Sb[M],[S[M][i]])
 ),
 if mod(i,3000)=0 then print("---",i,"/",length(S[M]),"is done")
);

/* now we obtain the result for n=M */
print("#S_A(",M,")=",length(Sa[M]),",#S_B(",M,")=",length(Sb[M]),
      ",#S(",M,")=",length(Sa[M])+length(Sb[M]),"=",
      binomial(M,ceiling(M/2)));

/* time to start induction */

for ct:M-1 thru 1 step -1 do(
 print("#S(",ct,")=",binomial(M,ct)*binomial(ct,ceiling(ct/2))),
/* S[ct]:MakeS(ct),*/
 S[ct]:MakeFiltS(fil,ct),
 turn:(-1)^ct, /* turn of which? */
 Sx:[],Sy:[],
 if turn=+1 then(Tx:Sa[ct+1]) else(Tx:Sb[ct+1]),

/*S0:[],
/*pairs:[[5,9],[9,5],[6,8],[8,6]],*/
/*pairs:[[6,8],[8,6]],*/ /* o wins */
/*pairs:[[6,16],[16,6],[7,15],[15,7],[8,14],[14,8],[10,12],[12,10]],*/
pairs:[[6,8],[8,6],[10,12],[12,10],[14,16],[16,14]],
sgn:-1,
for i:1 thru length(S[ct]) do(
 t0:S[ct][i],
 if turn=+1*sgn then(
  p:0,q:0,
  for j:1 thru length(pairs) do(
   if t0[pairs[j][1]]#0 and t0[pairs[j][1]]=t0[pairs[j][2]] then p:p+2,
   if t0[pairs[j][1]]=+1*sgn and t0[pairs[j][2]]=0 then p:p+1,
   if t0[pairs[j][1]]=-1*sgn and t0[pairs[j][2]]=0 then q:q+1
  ),
  if p<=1 and q<=1 then S0:append(S0,[t0])
 )else(
  p:0,
  for j:1 thru length(pairs) do(
   if t0[pairs[j][1]]#0 and t0[pairs[j][1]]=t0[pairs[j][2]] then p:p+2,
   if t0[pairs[j][1]]=+1*sgn and t0[pairs[j][2]]=0 then p:p+1,
   if t0[pairs[j][1]]=-1*sgn and t0[pairs[j][2]]=0 then p:p+2
  ),
  if p<=1 then S0:append(S0,[t0])
 )
),
S[ct]:S0,*/

 for i:1 thru length(S[ct]) do(
  flag:false,
  e:0,
  for j:1 thru M do(
   if S[ct][i][j]=0 then(
    t:copylist(S[ct][i]),
    t[j]:turn,
    if member(t,Tx) then(flag:true,j:M+1)
   )
  ),
  if flag=true then(
   Sx:append(Sx,[S[ct][i]])
  )else(
   Sy:append(Sy,[S[ct][i]])
  ),if mod(i,10000)=0 then print("---",i,"/",length(S[ct]),"is done")
 ),
 if turn=+1 then(Sa[ct]:Sx,Sb[ct]:Sy) else(Sb[ct]:Sx,Sa[ct]:Sy),
 print("#S_A(",ct,")=",length(Sa[ct]),",#S_B(",ct,")=",length(Sb[ct]),
       ",#S(",ct,")=",length(Sa[ct])+length(Sb[ct]),"=",
       binomial(M,ct)*binomial(ct,ceiling(ct/2)))
);
/*
for ct:M-1 thru 1 step -1 do(
 Sx:[],Sy:[],
 turn:(-1)^ct, /* turn of which? */
 if turn=+1 then(T:Sa[ct+1]) else(T:Sb[ct+1]),
 for i:1 thru length(T) do(
  for j:1 thru M do(
   if T[i][j]=turn then(
    t:copylist(T[i]),
    t[j]:0,
    if not member(t,Sx) then(Sx:append(Sx,[t]))
   )
  ) ,if mod(i,1000)=0 then print("Sa[",ct,"]",i,"/",length(T))
 ),
 if turn=-1 then(T:Sa[ct+1]) else(T:Sb[ct+1]),
 for i:1 thru length(T) do(
  for j:1 thru M do(
   if T[i][j]=turn then(
    t:copylist(T[i]),
    t[j]:0,
    if (not member(t,Sx)) and (not member(t,Sy)) then(Sy:append(Sy,[t]))
   )
  ) ,if mod(i,1000)=0 then print("Sb[",ct,"]",i,"/",length(T))
 ),
 if turn=+1 then(Sa[ct]:Sx,Sb[ct]:Sy) else(Sb[ct]:Sx,Sa[ct]:Sy),
 print("#S_A(",ct,")=",length(Sa[ct]),",#S_B(",ct,")=",length(Sb[ct]),
       ",#S(",ct,")=",length(Sa[ct])+length(Sb[ct]),"=",
       binomial(M,ct)*binomial(ct,ceiling(ct/2)))
);
*/
* * * * */

ExpCells(t):=block(
 [Sa0,Sb0,p,turn,i,t1,a,b],
 print(t),PrintTab(t),
 Sa0:Filt(t,Sa[M]),print(length(Sa0),"/",length(Sa[M])),
 Sb0:Filt(t,Sb[M]),print(length(Sb0),"/",length(Sb[M])),print("---"),
 p:makelist(0,i,1,M),
 turn:0,
 for i:1 thru M do(turn:turn+t[i]),
 turn:-2*turn+1,
 for i:1 thru M do(
  if t[i]=0 then(
   t1:makelist(0,i,1,M),
   t1[i]:turn,PrintTab(t+t1),
   a:length(Filt(t1,Sa0)),
   b:length(Filt(t1,Sb0)),
   p[i]:float(floor(float(a/(a+b))*1000)/1000),print(i,":",p[i]),print("---")
  )
 ),
 return(p)
);

/*********/
/* print */
/*********/

Stone(x):=block(
 if x=0 then(return("_"))
 else if x=1 then(return("o"))
 else(return("x"))
);
PrintTab(t):=block(
 [str,p,r],
 for i:1 thru N do(
  str:"",
  for j:1 thru N+1-i do(
   str:concat(str," ")
  ),
  for j:1 thru N do(
   str:concat(str,Stone(t[(i-1)*N+j])," ")
  ),
  if i=1 then(str:concat(str,Stone(t[M])," "))
  else(str:concat(str,Stone(t[(N-i+1)*N+1])," ")),
  print(str)
 ),
 str:concat(Stone(t[M])," "),
 for j:1 thru N do(
  str:concat(str,Stone(t[N+1-j])," ")
 ),
 print(str),
 print(""),
 return()
);
PrintTabs(T):=for i:1 thru length(T) do PrintTab(T[i]);

/*****************/
/* better place? */
/*****************/

AutoTurnStart(player):=block(
 Sa0:[],Sb0:[],
 ttt:makelist(0,i,1,M),
 if player=-1 then ttt[M]:1,
 print(ttt),PrintTab(ttt),
 print("you =",player)
);
AutoTurn(x,y):=block(
 [p,turn,t1,a,b,z,i,i0,r,xy,s],
 turn:0,
 for i:1 thru M do(turn:turn+ttt[i]),
 turn:-2*turn+1, /* turn of human */
 if not([x,y]=[0,0] or (1<=x and x<=N+1 and 1<=y and y<=N+1)) then(
  return("invalid input - out of board")
 ),
 if [x,y]=[0,0] or [x,y]=[1,N+1] or [x,y]=[N+1,1] then(
  z:M
 )elseif [x,y]=[N+1,N+1] then(
  z:1
 )else(
  xy:NormalizeXY([x-1,y-1]),x:xy[1]+1,y:xy[2]+1,
  z:N*(x-1)+y
 ),
 if ttt[z]#0 then return("invalid input - duplication"),
 ttt[z]:turn,
/* print(ttt),*/
 PrintTab(ttt),
/* if length(Sa0)+length(Sb0)=0 then(Sa0:Sa[M],Sb0:Sb[M]),
 Sa0:Filt(ttt,Sa0),
 Sb0:Filt(ttt,Sb0),
 print("#S =",length(Sa0),"+",length(Sb0)),print("---"),
 p:makelist(0,i,1,M),
 for i:1 thru M do(
  if ttt[i]=0 then(
   t1:makelist(0,i,1,M),
   t1[i]:-turn,
   a:length(Filt(t1,Sa0)),
   b:length(Filt(t1,Sb0)),
   p[i]:float(a/(a+b))
  )
 ),
 p:(-turn)*p,
 r:-2,
 for i:1 thru M do(
  if ttt[i]=0 and r<p[i] then(i0:i,r:p[i0])
 ),
 if r=-2 then return(),*/
 s:0,
 for i:1 thru M do(if ttt[i]=0 then s:s+1),
 if s>0 then(
  i0:MonteCarlo(ttt),
  ttt[i0]:(-turn),
/* print(ttt),*/
  print(""),PrintTab(ttt),
/* print(i0,p),*/
  print("COM : (",((i0-1)-mod(i0-1,N))/N+1,",",mod(i0-1,N)+1,")")
 ),
 print("Euler numbers =>",EulerNumber(ttt,+1),":",EulerNumber(ttt,-1)),
/* print("probability of win of 'o' =",abs(r)),*/
 return(i0)
);

/*                    */
/* Monte Carlo method */
/*                    */
R1:30;  /* basic */
R2:10;  /* top R2 places go to 2nd trial */
R3:300; /* 2nd */
MonteCarlo(t):=block(
 [turn,i,places,probs,tops,v,s,tab],
 turn:0,
 for i:1 thru M do turn:turn+t[i],
 turn:1-2*turn,
 places:[],
 for i:1 thru M do(
  if t[i]=0 then places:append(places,[i])
 ),
 probs:makelist(0,i,1,length(places)),
 for i:1 thru length(places) do(
  t[places[i]]:turn,
  probs[i]:ProbMC(t,-turn,R1),
  t[places[i]]:0
 ),
 tops:[],
 if length(places)<=R2 then(
  tops:makelist(0,i,1,length(places)),
  for i:1 thru length(tops) do(
   tops[i]:places[i]
  )
 )else(
  while length(tops)<R2 do(
   s:1,
   v:probs[1],
   for i:2 thru length(places) do(
    if (probs[i]-v)*turn>0 then(
     s:i,
     v:probs[i]
    )
   ),
   tops:append(tops,[places[s]]),
   probs[s]:-2*turn
  )
 ),
 tab:makelist(0,i,1,M),
 for i:1 thru length(tops) do tab[tops[i]]:+1,
 PrintTab(tab),
 probs:makelist(0,i,1,length(tops)),
 for i:1 thru length(tops) do(
  t[tops[i]]:turn,
  probs[i]:ProbMC(t,-turn,R3),
  t[tops[i]]:0
 ),
 print(tops,probs),
 s:1,
 v:probs[1],
 for i:2 thru length(tops) do(
  if (probs[i]-v)*turn>0 then(
   s:i,
   v:probs[i]
  )
 ),
 print("prob ~",v),
 return(tops[s])
);

/* computing probability (r times) */
ProbMC(t,turn,r):=block(
 [places,i,j,k,pp,q,tt,win],
 win:0,
 places:[],
 for i:1 thru M do(
  if t[i]=0 then places:append(places,[i])
 ),
 for k:1 thru r do(
  tt:copylist(t),
  pp:copylist(places),
  q:turn,
  while length(pp)>0 do(
   j:random(length(pp))+1,
   tt[pp[j]]:q,
   pp:delete(pp[j],pp),
   q:-q
  ),
  if EulerNumber(tt,+1)>0 then win:win+1
 ),
 return(float(win/r))
);