(* Glaesers Dominos                                                                  *)
(* Mathematica program by Oliver Knill, 31. December 2005                            *)

(* Fit the 28 Dominos {i,j},0 <=i<n,0<=j<n, n=7 onto the given matrix                *)
(* It takes about an hour to dig through all the possible cases                      *)
(* 128 = 2^7 solutions are found but they are all the same because of the orientation*)
(* of the seven (i,i) dominos                                                        *)
(* The program also verifies that there is a unique solution                         *)
(* and Glaesers matrix is indeed a Glaeser matrix (www.mathematik.com/Glaeser).      *)
(* I don't know whether there exist infinitely many Glaeser matrices                 *)

(* Algorithm: the program does what puzzle solvers do too. Put dominos until it      *)
(* is no more possible, then rearrange the last domino, try again etc                *)
(* Searching through this tree is done recursively. The procedure FindSolution(A,k)  *)
(* finds all possible places, where one can put a domino on the matrix A. Already    *)
(* occupied places of the matrix are labeled with the number 9                       *)
(* if all entries become 9's, the puzzle is solved.                                  *)

(* Usage: save this file glaesersdomino.m  and run it in Mathematica                 *)
(* It is reporting progress on the solutions, so that it is best run from a shell    *)
(* like    math< glaesersdomino.m     in linux or                                    *)
(* /Applications/Mathematica\ 5.1.app/Contents/MacOS/MathKernel <glaesersdomino      *)
(* in OSX. (Wolfram research does unfortunately not yet include an alias to the      *)
(* Mathkernel in OS X as in other Unix flavors)                                      *)

A={{3,6,2,0,0,4,4},
   {6,5,5,1,5,2,3},
   {6,1,1,5,0,6,3},
   {2,2,2,0,0,1,0},
   {2,1,1,4,3,5,5},
   {4,3,6,4,4,2,2},
   {4,5,0,5,3,3,4},
   {1,6,3,0,1,6,6}}
B={{0,1},{5,3},{3,2},{4,2},{3,1},{5,5},{0,0},{2,2},{5,0},{6,1},{5,1},{6,5},{5,2},{4,0},
   {4,3},{5,4},{6,4},{4,1},{1,1},{3,0},{6,6},{4,4},{6,0},{2,0},{6,3},{3,3},{2,1},{6,2}};
solution=Table[0,{Length[B]}]; m=Length[A]; n=Length[A[[1]]]; 
FitDomino[AA_,xx_]:=Module[{s,i,j},s={};
  Do[ Do[
       If[j<n && xx[[1]]==AA[[i,j]] && xx[[2]]==AA[[i,j+1]],s=Append[s,{{i,j},{i,j+1}}]]; 
       If[i<m && xx[[1]]==AA[[i,j]] && xx[[2]]==AA[[i+1,j]],s=Append[s,{{i,j},{i+1,j}}]]; 
       If[j<n && xx[[2]]==AA[[i,j]] && xx[[1]]==AA[[i,j+1]],s=Append[s,{{i,j+1},{i,j}}]]; 
       If[i<m && xx[[2]]==AA[[i,j]] && xx[[1]]==AA[[i+1,j]],s=Append[s,{{i+1,j},{i,j}}]],
    {j,1,n}],{i,1,m}]; Return[{xx,s}]];
Places[AA_,BB_]:=Module[{},PP={};Do[ss=FitDomino[AA,BB[[k]]];PP=Append[PP,ss],{k,1,Length[BB]}]; PP];
PlaceNumbers[AA_,BB_]:=Module[{},pps=Places[AA,BB];Table[Length[pps[[k,2]]],{k,Length[pps]}]];
pn=PlaceNumbers[A,B];
PN[AA_,BB_]:=Module[{mn=Length[BB]},pn=PlaceNumbers[AA,BB]; 10^40+Sum[10^k pn[[mn-k]],{k,0,mn-1}]]
RemoveDomino[AA_,qq_]:=Module[{AAA},AAA=AA; 
     AAA[[qq[[1,1]],qq[[1,2]]]]=9; AAA[[qq[[2,1]],qq[[2,2]]]]=9; Return[AAA]];
FindSolution[AA_,kk_]:=Module[{bb,pp,K,ll,AAA},AAA=AA;
   If[kk==Length[B],Write["solution.txt",solution]; 
    solution1=solution; Print[solution]; Print[MatrixForm[AAA]]];
   {bb,pp}=FitDomino[AA,B[[kk]]]; PP={bb,pp}; K=Length[pp];
   If[K>0 && Product[pn[[j]],{j,kk+1,Length[pn]}] >0,Do[
       solution[[kk]]=ll;
       qqq=pp[[ll]]; Print[PN[AA,B]]; Print[10^40+ll*10^(Length[B]-kk),"       ",count];
       AAA=RemoveDomino[AA,qqq]; FindSolution[AAA,kk+1],
     {ll,1,K}],count=count+1;  (*,else no solution found *) ]
]

Start:=Module[{}, count=0; Print["\n"]; FindSolution[A,1] ]; 

PlotSolution:=Module[{}, s={}; eps={0.05,0.05};
  solution0={3,1,2,2,1,1,3,1,1,4,1,1,1,1,2,1,1,1,1,2,1,1,1,1,1,1,1,1}; 
  nm=Length[solution0]; A1=A; 
  Do[{bb,pp}=FitDomino[A1,B[[kkkk]]]; PP={bb,pp}; qqq=pp[[solution0[[kkkk]]]]; 
    A2=RemoveDomino[A1,qqq]; A1=A2; Print[MatrixForm[A1]];Print["\n"]; 
    qqqq={ {qqq[[1,2]],nm-qqq[[1,1]]},{qqq[[2,2]],nm-qqq[[2,1]]}}; qqq=qqqq;  (* matrix -> descartes*)
    If[qqq[[1,1]]>qqq[[2,1]] || qqq[[1,2]]> qqq[[2,2]], qqqq={qqq[[2]],qqq[[1]]}; qqq=qqqq]; (*order*)
    s=Append[s,{RGBColor[0,0,0],Rectangle[qqq[[1]],qqq[[2]]+{1,1}]}];
    s=Append[s,{Hue[kkkk/nm],Rectangle[qqq[[1]]+eps,qqq[[2]]+{1,1}-eps]}],
  {kkkk,nm}]; A1=A; 
  Do[{bb,pp}=FitDomino[A1,B[[kkkk]]]; PP={bb,pp}; qqq=pp[[solution0[[kkkk]]]];
    A2=RemoveDomino[A1,qqq]; A1=A2;
    qqqq={ {qqq[[1,2]],nm-qqq[[1,1]]},{qqq[[2,2]],nm-qqq[[2,1]]}}; qqq=qqqq;  (* matrix -> descartes*)
    s=Append[s,{RGBColor[0,0,0],{Text[bb[[1]],qqq[[1]]+{1,1}/2],
                                 Text[bb[[2]],qqq[[2]]+{1,1}/2]}}],
  {kkkk,nm}]; S=Show[Graphics[s],AspectRatio->1,Frame->False,Axes->False]; 
  $TextStyle={FontFamily ->"Swiss721BT-Roman", FontSize->24};
  Display["glaesersdomino.ps",S,"EPS"];
  $TextStyle={FontFamily ->"Swiss721BT-Roman", FontSize->72};
  Display["glaesersdomino.gif",S,"GIF",ImageSize->1200]];

GifName[n_]:=Block[{u,v,w}, u=IntegerDigits[n]+48; v=ToCharacterCode["domino"];
  w=ToCharacterCode[".gif"]; FromCharacterCode[Join[v,u,w]]];

PlotDominos:=Module[{}, eps={0.05,0.05}; nm=Length[B]; 
  $TextStyle={FontFamily ->"Swiss721BT-Roman", FontSize->24};
  Do[ s={};  
   (* s={{RGBColor[0,0,0],Rectangle[{0,-1},{1,0}    ]}};
      s=Append[s,{RGBColor[0,0,0],Rectangle[{1,-1},{2,0}    ]}]; *)
      s=Append[s,{Hue[kkk/nm],    Rectangle[{0,-1}+eps,{1,0}-eps]}];
   (* s=Append[s,{Hue[kkk/nm],    Rectangle[{1,-1}+eps,{2,0}-eps]}];
      s=Append[s,{RGBColor[0,0,0],{Text[B[[kkk,1]],{0,-1}+{1,1}/2],
                                   Text[B[[kkk,2]],{1,-1}+{1,1}/2]}}]; *)
      S=Show[Graphics[s],AspectRatio->0.5,Frame->False,Axes->False,Background->Hue[kkk/nm]];
      Display[GifName[1000+kkk],S,"GIF", ImageSize->10], {kkk,1,nm}]; 
];


PlotProblem:=Module[{}, s={}; eps={0.05,0.05};
  s={RGBColor[0,0,0],Table[Text[A[[i,j]],{j,Length[A]-i}+{1,1}/2],{j,Length[A[[1]]]},{i,Length[A]}]}; 
  S=Show[Graphics[s],AspectRatio->1,Frame->False,Axes->False];
  $TextStyle={FontFamily ->"Swiss721BT-Roman", FontSize->24};
  Display["glaesersproblem.ps",S,"EPS"];
  $TextStyle={FontFamily ->"Swiss721BT-Roman", FontSize->72};
  Display["glaesersproblem.gif",S,"GIF",ImageSize->1200]];


Timing[Start]
(* 
PlotProblem
PlotDominos
PlotSolution
*)