(***************************************************
 *                                                 *
 *     MATHEMATICA PROGRAM BRINKHUISTRIPLES.M      *
 *                                                 *
 *     CONTENT:                                    *
 *     This program contains definitions to        *
 *     check generalised Brinkhuis triples.        *
 *     It accompanies the article "Improved        *
 *     bounds on the number of ternary             *
 *     square-free words" by Uwe Grimm.            *
 *     The Brinkhuis triples derived in the        *
 *     paper are explicitly contained.             *
 *                                                 *
 *     The program provides two main routines,     *
 *     BrinkhuisTripleCheck and                    * 
 *     BrinkhuisSpecialTripleCheck. The latter     *
 *     exploits the structure of a special         *
 *     Brinkhuis triple to reduce the number       *
 *     of words that have to be checked for        *
 *     square-freeness.                            *
 *                                                 *
 *     If the program is executed in the           *
 *     form, it starts to check the special        *
 *     Brinkhuis triples presented in the          *
 *     paper mentioned above. You may have         *
 *     to be patient if you are interested         *
 *     in the result for the largest triples,      *
 *     as it will take quite a while (several      *
 *     days on a Pentium 1GHz machine) to          *
 *     complete all checks.                        *
 *                                                 *
 *     DISCLAIMER:                                 *
 *     The program has been tested by the          *
 *     author, but no guarantee can be given       *
 *     that results are indeed correct - so        *
 *     use at your own risk!                       *
 *                                                 *
 *     AUTHOR:                                     *
 *     Uwe Grimm                                   *
 *     Applied Mathematics Department              *
 *     Faculty of Mathematics and Computing        *
 *     The Open University                         *
 *     Walton Hall                                 *
 *     Milton Keynes MK7 6AA                       *
 *     UK                                          *
 *                                                 *
 *     Please send any comments, bug reports,      *
 *     etc. to: u.g.grimm@open.ac.uk               *
 *                                                 *
 *     Version 1.0                                 *
 *     August 1, 2001                              *
 *                                                 *
 ***************************************************)

(*  DEFINITIONS  *)

Clear[ReverseWord,PermuteWord];
ReverseWord[w_String] := 
     StringReverse[w];
PermuteWord[w_String] := 
     StringReplace[w, {"0"->"1","1"->"2","2"->"0"}];

Clear[SubWords,NoSquare,SquareFree];
SubWords[w_String,len_Integer] := 
     Union[Table[StringTake[w,{i,i+len-1}],
                 {i,StringLength[w]-len+1}]];
NoSquare[w_String] /; Mod[StringLength[w],2]==0 := 
     (StringTake[#,StringLength[#]/2]=!=
      StringTake[#,-StringLength[#]/2])&[w];
SquareFree[w_String] := 
     Apply[And,Table[Apply[And,Map[NoSquare,SubWords[w,k]]],
                     {k,2,StringLength[w],2}]];

Clear[BrinkhuisSpecialTriple];
BrinkhuisSpecialTriple[www_List] :=
NestList[Map[PermuteWord,#]&,Union[www,Map[ReverseWord,www]],2];

Clear[BrinkhuisTripleCheck];
BrinkhuisTripleCheck[www_List] :=
   Module[{threewords={"010","012","020","021",
                       "101","102","120","121",
                       "201","202","210","212"},
           ww,wl1,wl2,wl3,lw1,lw2,lw3,
           i,j,k,n,nn,tn=0,sqfree=True,sqf},
          If[Length[Union[Flatten[www]]]=!=Length[Flatten[www]],
             Print["triple contains identical words"];
             Return[False]];	  
          If[Length[Union[Map[StringLength,Flatten[www]]]]=!=1,
             Print["triple contains words of different lengths"];
             Return[False]];      
          Print["checking a triple consisting of (",
                Length[www[[1]]],",",Length[www[[2]]],",",
                Length[www[[3]]],") words of length ",
                Union[Map[StringLength,Flatten[www]]][[1]]]; 
          Do[ww=threewords[[n]];
             wl1 = www[[1+ToExpression[StringTake[ww,{1,1}]]]];
             wl2 = www[[1+ToExpression[StringTake[ww,{2,2}]]]];
             wl3 = www[[1+ToExpression[StringTake[ww,{3,3}]]]];
             lw1 = Length[wl1];
             lw2 = Length[wl2];
             lw3 = Length[wl3];
             nn  = lw1*lw2*lw3;
             tn += nn;
             sqf = Apply[And,Map[SquareFree,
                         Flatten[Table[StringJoin[wl1[[i]],
                                                  wl2[[j]],
                                                  wl3[[k]]],
	 	                       {i,lw1},
                                       {j,lw2},
                                       {k,lw3}]]]];
             Print[" - checked ",lw1*lw2*lw3," words of type ",
                   ww,": ",sqf];
             sqfree = And[sqfree,sqf],
	    {n,Length[threewords]}];
          Print["all checks completed, a total of ",tn,
                " words were checked"];
	  If[sqfree,
	     Print["no squares found - this is indeed a ",
                   "Brinkhuis triple!"],
             Print["squares detected - this is not a ",
                   "Brinkhuis triple!"]];
          sqfree];

Clear[BrinkhuisSpecialTripleCheck];
BrinkhuisSpecialTripleCheck[wws_List] :=
   Module[{threewords={"010","012","020"},
           ww,wl1,wl2,wl3,lw1,lw2,lw3,
           i,j,k,m,n,nn,tn=0,sqfree=True,sqf,diffword,
           www=BrinkhuisSpecialTriple[wws],
           npal=0,nword=Length[wws]},
	  npal=Count[Map[ReverseWord[#]===#&,wws],True];
          If[Length[Union[Flatten[www]]]=!=Length[Flatten[www]],
             Print["triple contains identical words"];
             Return[False]];	  
          If[Length[Union[Map[StringLength,Flatten[www]]]]=!=1,
             Print["triple contains words of different lengths"];
             Return[False]];      
          Print["checking a special triple consisting of (",
                Length[www[[1]]],",",Length[www[[2]]],",",
                Length[www[[3]]],") words of length ",
                Union[Map[StringLength,Flatten[www]]][[1]]];
          Print["signature of triple is (",
                npal,",",nword-npal,")"];
          Do[ww=threewords[[n]];
             wl1 = www[[1+ToExpression[StringTake[ww,{1,1}]]]];
             wl2 = www[[1+ToExpression[StringTake[ww,{2,2}]]]];
             wl3 = www[[1+ToExpression[StringTake[ww,{3,3}]]]];
             lw1 = Length[wl1];
             lw2 = Length[wl2];
             lw3 = Length[wl3];
             nn  = lw1*lw2*lw3;
             sqf = Union[Flatten[Table[StringJoin[wl1[[i]],
                                                  wl2[[j]],
                                                  wl3[[k]]],
	 	                       {i,lw1},
                                       {j,lw2},
                                       {k,lw3}]]];
             If[Length[sqf]=!=nn,
	        Print["Error detected - not all words replacing ",
                      ww," are different from each other"];];
             (* we don't need to check words that are mirror images *)
             Do[If[ReverseWord[sqf[[m]]]=!=sqf[[m]],
                   If[MemberQ[sqf,ReverseWord[sqf[[m]]]],
                      sqf=Delete[sqf,m]]],
                {m,Length[sqf],2,-1}];
             diffword = Length[sqf];
             tn += diffword;
             sqf = Apply[And,Map[SquareFree,sqf]];
             Print[" - checked ",diffword,
                   " non-equivalent words among ",
                   nn," words of type ",ww,": ",sqf];
             sqfree = And[sqfree,sqf],
	    {n,Length[threewords]}];
          Print["all checks completed, a total of ",tn,
                " words were checked"];
	  If[sqfree,
	     Print["no squares found - this is indeed a ",
                   "Brinkhuis triple!"],
             Print["squares detected - this is not a ",
                   "Brinkhuis triple!"]];
          sqfree];

(*  CHECKS OF SPECIAL BRINKHUIS TRIPLES PRESENTED IN THE PAPER  *)

(*  n=13  *)
Clear[specialbrinktrip13,brinktrip13];
specialbrinktrip13 = {"0121021201210"};
brinktrip13        = BrinkhuisSpecialTriple[specialbrinktrip13];
Print["checking special Brinkhuis 1-triple of length 13:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip13];
Print[" "];

(*  n=18  *)
Clear[specialbrinktrip18,brinktrip18];
specialbrinktrip18 = {"012021020102120210"};
brinktrip18        = BrinkhuisSpecialTriple[specialbrinktrip18];
Print["checking special Brinkhuis 2-triple of length 18:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip18];
Print[" "];

(*  n=23  *)
Clear[specialbrinktrip23,brinktrip23];
specialbrinktrip23 = {"01210201021012021201210",
                      "01210212021012021201210"};
brinktrip23        = BrinkhuisSpecialTriple[specialbrinktrip23];
Print["checking special Brinkhuis 3-triple of length 23:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip23];
Print[" "];

(*  n=25  *)
Clear[specialbrinktrip25,brinktrip25];
specialbrinktrip25 = {"0121020102101201021201210",
                      "0121021201021012021201210",
		      "0121021202102012021201210"};
brinktrip25        = BrinkhuisSpecialTriple[specialbrinktrip25];
Print["checking special Brinkhuis 5-triple of length 25:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip25];
Print[" "];

(*  n=29, case 1  *)
Clear[specialbrinktrip29a,brinktrip29a];
specialbrinktrip29a = {"01202102012101202120102120210",
                       "01202120102012021012102120210",
                       "01202120102012021020102120210",
                       "01202120121012021012102120210"};
brinktrip29a        = BrinkhuisSpecialTriple[specialbrinktrip29a];
Print["checking special Brinkhuis 6-triple of length 29:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip29a];
Print[" "];

(*  n=29, case 2  *)
Clear[specialbrinktrip29b,brinktrip29b];
specialbrinktrip29b = {"01210201021201020121021201210",
                       "01210201021202101201021201210",
                       "01210201021202102012021201210"};
brinktrip29b        = BrinkhuisSpecialTriple[specialbrinktrip29b];
Print["checking special Brinkhuis 6-triple of length 29:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip29b];
Print[" "];

(*  n=30  *)
Clear[specialbrinktrip30,brinktrip30];
specialbrinktrip30 = {"012102010210120102012021201210",
                      "012102010212012102012021201210",
                      "012102010212021020121021201210",
                      "012102120210120102012021201210"};
brinktrip30        = BrinkhuisSpecialTriple[specialbrinktrip30];
Print["checking special Brinkhuis 8-triple of length 30:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip30];
Print[" "];


(*  n=33  *)
Clear[specialbrinktrip33,brinktrip33];
specialbrinktrip33 = {"012021020121012010212012102120210",
                      "012021020121021201021012102120210",
                      "012021020121021201210120102120210",
                      "012021201020120210121020102120210",
                      "012021201020121012021012102120210",
                      "012021201021012010212012102120210"};
brinktrip33        = BrinkhuisSpecialTriple[specialbrinktrip33];
Print["checking special Brinkhuis 12-triple of length 33:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip33];
Print[" "];

(*  n=35  *)
Clear[specialbrinktrip35,brinktrip35];
specialbrinktrip35 = {"01202102010212010201202120102120210",
                      "01202102010212010201210120102120210",
                      "01202102012101201020121020102120210",
                      "01202102012101202120121020102120210",
                      "01202102012102120210121020102120210",
                      "01202120102012101201021012102120210",
                      "01202120102012102010210120102120210",
                      "01202120102120210201021012102120210",
                      "01202120102012102120121020102120210",
                      "01202120102101210201210120102120210"};
brinktrip35    = BrinkhuisSpecialTriple[specialbrinktrip35];
Print["checking special Brinkhuis 18-triple of length 35:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip35];
Print[" "];

(*  n=36  *)
Clear[specialbrinktrip36,brinktrip36];
specialbrinktrip36 = {"012102010210120212010210121021201210",
                      "012102010210120212012101201021201210",
                      "012102010210121021201020121021201210",
                      "012102010210121021201021012021201210",
                      "012102010210121021202101201021201210",
                      "012102010212012101201020121021201210",
                      "012102010212012101201021012021201210",
                      "012102010212012102120210121021201210",
                      "012102010212021012010210121021201210",
                      "012102010212021020102101201021201210",
                      "012102120102101202120102012021201210",
                      "012102120102101210201021012021201210",
                      "012102120102101210212021012021201210",
                      "012102120121012010212021012021201210",
                      "012102120121012021201021012021201210",
		      "012102120121020102120102012021201210"};
brinktrip36        = BrinkhuisSpecialTriple[specialbrinktrip36];
Print["checking special Brinkhuis 32-triple of length 36:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip36];
Print[" "];

(*  n=40  *)
Clear[specialbrinktrip40,brinktrip40];
specialbrinktrip40 = {"0120210201210120102120210121020102120210",
                      "0120210201210120212010210121020102120210",
                      "0120210201210212010210120212012102120210",
                      "0120210201210212012101201021012102120210",
                      "0120210201210212012101202120121020120210",
                      "0120210201210212012102010210120102120210",
                      "0120210201210212012102010212012102120210",
                      "0120210201210212012102012021012102120210",
                      "0120210201210212012102012021020102120210",
                      "0120210201210212021020120212012102120210",
                      "0120212010201202101210201021012102120210",
                      "0120212010201210120102012021012102120210",
                      "0120212010201210120102101202120102120210",
                      "0120212010201210120102120121020102120210",
                      "0120212010201210120210201021012102120210",
                      "0120212010201210120210201202120102120210",
                      "0120212010201210120212010210120102120210",
                      "0120212010201210212010201202120102120210",
                      "0120212010201210212012101202120102120210",
                      "0120212010210120102012101202120102120210",
                      "0120212010212021012021201021012102120210",
                      "0120212010212021012102010212012102120210",
                      "0120212010212021012102120102012102120210",
                      "0120212010212021020102120102012102120210"};
brinktrip40        = BrinkhuisSpecialTriple[specialbrinktrip40];
Print["checking special Brinkhuis 48-triple of length 40:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip40];
Print[" "];

(*  n=41  *)
Clear[specialbrinktrip41,brinktrip41];
specialbrinktrip41 = {"01202102012101201021202101202120102120210",
                      "01202102012101202120102012021012102120210",
                      "01202102012101202120102012021020102120210",
                      "01202102012102120102012101202120102120210",
                      "01202102012102120121012010212012102120210",
                      "01202102012102120121020120212012102120210",
                      "01202120102012021012010210121020102120210",
                      "01202120102012021012102010210120102120210",
                      "01202120102012021012102010212012102120210",
                      "01202120102012021012102012021020102120210",
                      "01202120102012021012102120102012102120210",
                      "01202120102012021012102120121020102120210",
                      "01202120102012021020102120102012102120210",
                      "01202120102012021020102120121020102120210",
                      "01202120102012101201020120212012102120210",
                      "01202120102012101202102010210120102120210",
                      "01202120102012101202102010212012102120210",
                      "01202120102012101202102012021012102120210",
                      "01202120102012102120102012021012102120210",
                      "01202120102012102120102101202120102120210",
                      "01202120102012102120121012021012102120210",
                      "01202120102012102120210201021012102120210",
                      "01202120102012102120210201202120102120210",
                      "01202120102101201020121012021012102120210",
                      "01202120102101201021202101202120102120210",
                      "01202120102120210121021201021012102120210",
                      "01202120102120210201202120102012102120210",
                      "01202120121012010201202120102012102120210",
                      "01202120121012021012102010212012102120210",
                      "01202120121012021012102120102012102120210",
                      "01202120121012021020102120102012102120210",
                      "01202102012102120210201202120121020120210",
                      "01202120121012010201210201021012102120210",
                      "01202120121021201021012010212012102120210"};
brinktrip41        = BrinkhuisSpecialTriple[specialbrinktrip41];
Print["checking special Brinkhuis 65-triple of length 41:"];
BrinkhuisSpecialTripleCheck[specialbrinktrip41];
Print[" "];

(*  END OF PROGRAM BRINKHUISTRIPLES.M  *)

Quit[];
