##########################################################
#
#
#
#                merlin.mpl   simulates lights out and solves it
#
#
#
#wdj,5-2000, 12-2002: bug in play_keychain found by Jaap S.
###########################################################

elemats:=proc(n)
local i,j:
global E;
for i from 1 to n do
for j from 1 to n do
E[i,j,n]:=array(sparse,1..n,1..n):
E[i,j,n][i,j]:=1:
od:
od:
end:

convert2list:=proc(A::array)
local m,n,i,j,v,A0;
m:=rowdim(A);
v:=[];
A0:=convert(A,listlist):
for i from 1 to m do
 v:=[op(v),op(A0[i])];
od;
RETURN(v):
end:


convert2array:=proc(v::list,m)
local N,i,j,A0;
N:=nops(v);
A0:=array(sparse,1..m,1..N/m):
for i from 1 to m do
 for j from 1 to N/m do
  A0[i,j]:=v[j+(i-1)*N/m];
 od;
od;
RETURN(convert(A0,matrix)):
end:

solve5x5:=proc(S::array)
#S is the 5x5 state matrix
#not all such matrices are solvable
local a,b,i,j,v,B,I4,O4,A5,sol,soln,sol0:
B:=BandMatrix([1,1,1], 1, 5, outputoptions=[storage=rectangular]);
I4:=BandMatrix([0,1,0], 1, 5, outputoptions=[storage=rectangular]);
O4:=BandMatrix([0,0,0], 1, 5, outputoptions=[storage=rectangular]);
A5:=blockmatrix(5,5,[B,I4,O4,O4,O4,I4,B,I4,O4,O4,O4,I4,B,I4,O4,O4,O4,I4,B,I4,O4,O4,O4,I4,B]);
v:=convert2list(S);
sol:=Linsolve(evalm(A5),convert(v,array)) mod 2;
#sol:=linsolve(convert(A5,matrix),convert(v,vector));
if sol=NULL then print(`No solution for `,S); RETURN(`try again`); fi;
sol0:=(a,b)->[seq(subs({_t[24]=a,_t[25]=b},sol[j]),j=1..25)]:
#print(11,sol0(a,b),sol);
soln:=[seq(sol0(1,0)[i] mod 2,i=1..25)];
#print(22,soln);
RETURN(convert2array(soln,5));
end:

play_5x5:=proc(S::array,move::list)
#move=[i,j] is the move in matrix notation, so [1,1] is the
#top left corner, [5,1] is the bottom left corner, ...
#s is the state array
local a,b,i,j,v,B,I4,O4,A5,Snew,sol,soln,sol0:
B:=BandMatrix([1,1,1], 1, 5, outputoptions=[storage=rectangular]);
I4:=BandMatrix([0,1,0], 1, 5, outputoptions=[storage=rectangular]);
O4:=BandMatrix([0,0,0], 1, 5, outputoptions=[storage=rectangular]);
A5:=blockmatrix(5,5,[B,I4,O4,O4,O4,I4,B,I4,O4,O4,O4,I4,B,I4,O4,O4,O4,I4,B,I4,O4,O4,O4,I4,B]);
v:=convert2list(S);
Snew:=convert2list(evalm(S+convert2array(convert(row(A5,(move[1]-1)*5+move[2]),list),5))) mod 2:
Snew:=convert2array(Snew,5):
RETURN(evalm(Snew)):
end:

solve6x6:=proc(S::array)
#S is the 6x6 state matrix
local i,j,v,B,I4,O4,A6,sol,Snew,soln:
B:=BandMatrix([1,1,1], 1, 6, outputoptions=[storage=rectangular]);
I4:=BandMatrix([0,1,0], 1, 6, outputoptions=[storage=rectangular]);
O4:=BandMatrix([0,0,0], 1, 6, outputoptions=[storage=rectangular]);
A6:=blockmatrix(6,6,[B,I4,O4,O4,O4,O4,I4,B,I4,O4,O4,O4,O4,I4,B,I4,O4,O4,O4,O4,I4,B,I4,O4,O4,O4,O4,I4,B,I4,O4,O4,O4,O4,I4,B]);
v:=convert2list(S);
sol:=linsolve(convert(A6,matrix),v);
#print(11,sol);
soln:=[seq(sol[i] mod 2,i=1..36)];
#print(22,soln);
RETURN(convert2array(soln,6));
end:

#the following procedure is, so far, not working
#for some unknown reason
solveorbix:=proc(v::list)
#S is the 12x12 state vector
local A,i,j,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,sol,soln:
B1:=[0,1,1,1,1,1,0,0,0,0,0,0]:
B2:=[1,0,1,0,0,1,1,0,0,0,1,0]:
B3:=[1,1,0,1,0,0,1,1,0,0,0,0]:
B4:=[1,0,1,0,1,0,0,1,1,0,0,0]:
B5:=[1,0,0,1,0,1,0,0,1,1,0,0]:
B6:=[1,1,0,0,1,0,0,0,0,1,1,0]:
B7:=[0,1,1,0,0,0,0,1,0,0,1,1]:
B8:=[0,0,1,1,0,0,1,0,1,0,0,1]:
B9:=[0,0,0,1,1,0,0,1,0,1,0,1]:
B10:=[0,0,0,0,1,1,0,0,1,0,1,1]:
B11:=[0,1,0,0,0,1,1,0,0,1,0,1]:
B12:=[0,0,0,0,0,0,1,1,1,1,1,0]:
A:=matrix([B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12]);
#print(A);
#latex(A);
sol:=linsolve(A,v);
if sol=NULL then print(`No solution for `,v); RETURN(`try again`); fi;
#print(11,sol);
soln:=[seq(sol[i] mod 2,i=1..12)];
#print(22,soln);
RETURN(soln);
end:

#returns codim of mxn merlin's machine
#soln space in state space over F_2
codim_merlin:=proc(m::integer,n::integer)
local rank_2,A,i,j,GA,k,l,move,mv;
for i from 1 to m do
 for j from 1 to n do
  move[i,j]:=array(sparse,1..m,1..n):
  for k from 1 to m do
   for l from 1 to n do
    if abs(i-k)+abs(j-l)<=1 then move[i,j][k,l]:=1: fi;
   od;
 od;
 mv[i,j]:=convert2list(move[i,j]):
 #print(i,j,evalm(move[i,j]),mv[i,j]);
 od;
od;
A:=matrix([seq(seq(mv[i,j],j=1..n),i=1..m)]);
GA:=Gausselim(A) mod 2;
rank_2:=0:
for i from 1 to rowdim(A) do
 if iszero(row(GA,i)) then rank_2:=rank_2+1; fi;
od:
RETURN(rank_2);
end:

#returns toggle matrix of mxn merlin's machine
#soln space in state space over F_2
togglemat_merlin:=proc(m::integer,n::integer)
local rank_2,A,i,j,GA,k,l,move,mv;
for i from 1 to m do
 for j from 1 to n do
  move[i,j]:=array(sparse,1..m,1..n):
  for k from 1 to m do
   for l from 1 to n do
    if abs(i-k)+abs(j-l)<=1 then move[i,j][k,l]:=1: fi;
   od;
 od;
 mv[i,j]:=convert2list(move[i,j]):
# print(i,j,evalm(move[i,j]),mv[i,j]);
 od;
od;
A:=matrix([seq(seq(mv[i,j],j=1..n),i=1..m)]);
RETURN(evalm(A));
end:

#solve the keychain lights out, s is a 4x4 state vector
solve_keychain:=proc(s::array)
local keychain_togglematrix,A,v0,i,j,k,sol,soln,B,I4,O4,A4;
elemats(16):
B:=BandMatrix([1,1,1], 1, 4, outputoptions=[storage=rectangular]);
I4:=BandMatrix([0,1,0], 1, 4, outputoptions=[storage=rectangular]);
O4:=BandMatrix([0,0,0], 1, 4, outputoptions=[storage=rectangular]);
A4:=blockmatrix(4,4,[B,I4,O4,O4,I4,B,I4,O4,O4,I4,B,I4,O4,O4,I4,B]);
keychain_togglematrix:=
evalm(A4+E[1,4,16]+E[1,13,16]+E[2,14,16]+E[3,15,16]+E[4,1,16]+E[4,16,16]+E[5,8,16]+E[8,5,16]+E[9,12,16]+E[12,9,16]+E[13,1,16]+E[13,16,16]+E[14,6,16]+E[15,3,16]+E[16,4,16]+E[16,13,16]);
 A:=keychain_togglematrix:
#latex(A);
 v0:=convert2list(s):
 soln:=Linsolve(A,convert(v0,array)) mod 2:
 if nops(soln)=1 then 
 sol:=[seq(subs({_t[15]=0,_t[16]=0},soln[k]),k=1..16)]:
 print(`state: `,evalm(s)):
 print(` soln: `):
 RETURN(convert2array(sol,4)):
fi:
RETURN(`Sorry, no solution. Try again.`):
end:

##This has a bug - in plays [4,3] incorrectly.
play_keychain:=proc(S::array,move::list)
#move=[i,j] is the move in matrix notation, so [1,1] is the
#top left corner, [5,1] is the bottom left corner, ...
#s is the state array
local a,b,i,j,v,B,A4,I4,O4,A,Snew,sol,soln,sol0:
elemats(16):
B:=BandMatrix([1,1,1], 1, 4, outputoptions=[storage=rectangular]);
I4:=BandMatrix([0,1,0], 1, 4, outputoptions=[storage=rectangular]);
O4:=BandMatrix([0,0,0], 1, 4, outputoptions=[storage=rectangular]);
A4:=blockmatrix(4,4,[B,I4,O4,O4,I4,B,I4,O4,O4,I4,B,I4,O4,O4,I4,B]);
A:=evalm(A4+E[1,4,16]+E[1,13,16]+E[2,14,16]+E[3,15,16]+E[4,1,16]+E[4,16,16]+E[5,8,16]+E[8,5,16]+E[9,12,16]+
E[12,9,16]+E[13,1,16]+E[13,16,16]+E[14,2,16]+E[15,3,16]+E[16,4,16]+E[16,13,16]);
#print(A);
v:=convert2list(S);
Snew:=convert2list(evalm(S+convert2array(convert(row(A,(move[1]-1)*4+move[2]),list),4))) mod 2:
Snew:=convert2array(Snew,4):
RETURN(evalm(Snew)):
end: