# What is the fastest way to count square-free words?

## Background

A word is a string of letters in an alphabet. A square-free word has no adjacent repeating substring. For example, (in the ternary alphabet of {0,1,2}) the words 00, 012121, and 0212012021 are all not square-free because they contain repeated subwords, but a word like 1201021 or 0102012021 is square-free.

Here is code to test if a string over the alphabet {0,1,2} is square-free, but it is too slow:

ReverseWord[w_String] :=
StringReverse[w];
PermuteWord[w_String] :=
StringReplace[w, {"0" -> "1", "1" -> "2", "2" -> "0"}];
SubWords[w_String, len_Integer] :=
Module[{i}, 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}]];

SquareFree["0101"] (* returns False *)
SquareFree["0102012021"] (* returns True *)


I'm sure there are some wizards out there that knows the fastest way to do this! Here is a big squarefree word to test your fast functions on:

benchmark = Uncompress@"1:eJzt0sEJgDAQRNGUoh3sbinagQ3Y/\
03ISYIQQoSvzJ5UIknm7azHue1lKcXcws2jPuP2bvUrHtdj8P92vd3n7fN7+3/9/rPnz/\
rQ9//7/Gl/Oj/dH9qf9sv5a/vT+en+0P60X85f25/OT/dH3Z/OT/dH3Z/OT/\
eH9qf9cv7a/nR+uj+0P+2X89f2p/PT/aH9ab+cv7Y/nZ/uj7o/nZ/uj7r/\
wP0uYjwiEQ==";


## Problem

The ultimate way to determine square-freeness of words is to enumerate them! So the real problem is to generate all ternary squarefree words of length n efficiently.

Here was my original implementation: the function ExtendSquareFreeQ[w] checks if word w is squarefree assuming that StringTake[w, 1;;-2] is squarefree. Then, starting from the list AllSquareFree4, we build up the next AllSquareFree5, and so forth...

iExtendSquareFreeQ = Compile[{{str, _Integer, 1}},
Module[{sl = Length[str]/2 // Floor, i},
For[i = 2, i <= sl, i++,
If[Take[str, {-2 i, -(i + 1)}] === Take[str, {-i, -1}],
Return[False]]
];
True
], CompilationTarget -> "C", RuntimeOptions -> "Speed",
RuntimeAttributes -> {Listable}, Parallelization -> True
];
ExtendSquareFreeQ[str_String] := iExtendSquareFreeQ[ToCharacterCode[str]];
DistributeDefinitions[ExtendSquareFreeQ];

ExtendAllSquareFree[list_] := Module[
{cand, sfcand},
cand = Flatten @ Map[
Function[w,
Function[s, StringJoin[w, s]] /@
Complement[{"0", "1", "2"}, {StringTake[w, -1]}]
],
list
];
sfcand = ParallelMap[ExtendSquareFreeQ, cand];
Pick[cand, sfcand]
];

AllSquareFree2 = {"01"};
temp = AllSquareFree2;
byteCounts = {};
times = {};
Monitor[Table[
{time, temp} = AbsoluteTiming[ExtendAllSquareFree[temp]];
temps = Compress@temp;
AppendTo[times, time];
AppendTo[byteCounts, ByteCount[temps]/1000000];, {n, 2, 111}],
ListPlot[{times, byteCounts}, Mesh -> Full,
PlotRange -> {Automatic, All}, Filling -> Axis,
AxesLabel -> {"n", None}, InterpolationOrder -> 2, Joined -> True,
PlotLabel -> Text@Row[{
"n: " <> ToString[n],
"count: " <> ToString[Length[temp]*6],
"time: " <> ToString[N[time, 2]] <> " s",
"mem: " <> ToString[N[byteCounts // Last, 2]] <> " mB"},
" | "]]]


But both the time and memory complexities are way too high for $n>60$.

Is there a much more efficient way to perform this search? Perhaps storing the lists of strings in special data structures and traverse them in fancy ways?

I'll take any solutions of linear or quadratic complexity, or that can reach $n >= 111.$

The main difficulty is that when storing all the words of length 63 Mathematica runs out of memory! I was hoping for some brilliant optimized data structures, but I guess I might just have to use a tree...

-
You can assume the first two elements are 01. Any sq-free strings starting that way will, on permutations of {1,2,3}, give rise to all others. Might save a factor of 6 in an all-out enumeration. –  Daniel Lichtblau Aug 22 '12 at 22:09
@DanielLichtblau I just added your suggestion :) –  M.R. Aug 22 '12 at 22:20
This gives a bound between 197 and 1.9*10^16 ternary squarefree words of length 111. I honestly hope the result is nearer the lower bound –  Rojo Aug 22 '12 at 22:39
@Rojo, this shows that the number will be something around 10^13. –  Simon Woods Aug 23 '12 at 8:40
Then, using 2 bits for every solution we would just need 2.3TB. If the list of strings is required, who knows, about 200TB –  Rojo Aug 23 '12 at 10:26

This has [quadratic, he said] actually perhaps cubic complexity in a worst case (okay, now I'm just confused. More below).. Not the fastest of the lot, but it seems reasonable, or at least not entirely unreasonable. Requires some thought for me to see what I'm doing that keeps it relatively slow.

squareFree[wrd_String] := squareFreeC[ToCharacterCode[wrd]]

notSqrFreeAtTail = Compile[{{chars, _Integer, 1}}, Module[
{j, k, n = Length[chars], nhalf},
nhalf = Floor[n/2];
For[j = 1, j <= nhalf, j++,
For[k = 0, k < j, k++,
If[chars[[n - k]] != chars[[n - j - k]], Break[]];
];
If[k == j, Return[True]];
];
False
]];

squareFreeC = Compile[{{chars, _Integer, 1}}, Module[
{j, n = Length[chars]},
For[j = 2, j <= n, j++,
If[notSqrFreeAtTail[chars[[1 ;; j]]], Return[False]];
];
True
], CompilationOptions -> {"InlineExternalDefinitions" -> True},
RuntimeOptions -> "Speed", CompilationTarget -> "C"];

Timing[squareFree[benchmark]]

(* Out[12]= {1.22, True} *)


Re complexity: If I did this correctly then the total number of character comparisons should be no more than n choose 2 with n the length of the input. Ergo quadratic complexity. With particular emphasis in that "if".

--- edit ---

This next code enumerates all square free ternary words up to length n. It explicitly maintains the word list (as lists of integers range 0-2). It also records sufficient information that one could reconstruct all smaller such words, and gives a count of the full set. I did not try for the absolute maximal compression for the method of recording new words, but I believe that can be done in such a way as to take up only a few bits per new square free word. The basic idea is, recursively, to store only the new last letter from a given word of length one less. What I do is store a value between 0 and 7 (actually 6). The binary bits that are set indicate the corresponding letters that can follow this particular word with the result still being square free.

--- edit #2 (changed code) ---

squareFreeTernaryWords[n_Integer] :=
Module[{jwds, allwds, count = 1, jwd, newjwds, jnew, newwd},
allwds = Table[{}, {n}];
allwds[[1]] = {0};
allwds[[2]] = {1};
jwds = {{0, 1}};
Do[
allwds[[j]] = ConstantArray[0, Length[jwds]];
newjwds = Reap[
Do[
jwd = jwds[[k]];
jnew = 0;
For[i = 0, i <= 2, i++,
If[i != Last[jwd],
newwd = Append[jwd, i];
If[squareFreeC[newwd],
count++;
jnew = jnew + 2^i;
Sow[newwd]];
]];
allwds[[j, k]] = jnew;
, {k, Length[jwds]}];][[2, 1]];
jwds = newjwds;
Print[{j, Length[jwds]}];,
{j, 3, n}];
{jwds, count}]


--- end edit #2 ---

I did not try to use Compile but it could be rewritten to do so. It will use whatever square-free tester you like, provided that it operates correctly when given a list of integers. It should be posible to parallelize the k loop.

--- end edit ---

--- Mike's edit ---

## Parallelization

I tried boiling down your code to just get the counts in parallel (for the k loop) for my 24 core machine. Replacing the ParallelDo[] with a Do[] makes it work, not sure why it isn't working...

SetSharedVariable[jwds, newjwds, count, jwd];
DistributeDefinitions[squareFree,notSqrFreeAtTail,squareFreeC];
SetSharedVariable[newjwds, count];
count = 1;
jwds = {{0, 1}};
Do[
newjwds = {};
ParallelDo[
jwd = jwds[[k]];
For[i = 0, i <= 2, i++,
If[i == Last[jwd], Continue[], newwd = Append[jwd, i]];
If[squareFreeC[newwd], count++; AppendTo[newjwds, newwd];];
];
,{k, Length[jwds]}
];

jwds = newjwds;
Print[{j, 6*Length[jwds]}];
,{j, 3, 111}
]


--- end Mike's edit ---

--- edit #3 ---

I confess I had no joy from ParallelXXX. Everything I tried either did not work or was hugely slower than a straight serial computation.

If I jettison the tree emulation (and thus give up on recovering all words generated) then I can readily get the full body of code through Compile. I show one alternative below. It uses Leonid's square free tester because that was the fastest I saw in this thread. I do not think all the option settings are needed but kept them just in case.

With[{squareFreeQLSC = squareFreeQLSC},
squareFreeTernaryWordsC =
Compile[{{n, _Integer}},
Module[{jwds, count = 1, newjwds, jwd, newwd, children, i, m},
jwds = {{0, 1}};
Do[
newjwds = Table[
jwd = jwds[[k]];
children = ConstantArray[3, 2 j];
m = 0;
For[i = 0, i <= 2, i++,
If[i != Last[jwd],
newwd = Append[jwd, i];
If[squareFreeQLSC[newwd]==0,
(*children[[m+1;;m+j]]=newwd;m+=j*)

Do[children[[++m]] = newwd[[l]], {l, j}]];
];
];
children
, {k, Length[jwds]}];
jwds = Partition[Flatten[newjwds], j];
jwds = Select[jwds, FreeQ[#, 3] &];
count += Length[jwds];
Print[{j, count}];
,
{j, 3, n}];
count
], CompilationOptions -> {"InlineExternalDefinitions" -> True,
"InlineCompiledFunctions" -> True}, CompilationTarget -> "C",
RuntimeOptions -> "Speed"]];


This next takes up substantial memory. But it is fairly fast.

In[6]:= AbsoluteTiming[sf = squareFreeTernaryWordsC[55];]

(* During evaluation of In[6]:= {3,3}

During evaluation of In[6]:= {4,6}

During evaluation of In[6]:= {5,11}

During evaluation of In[6]:= {6,18}

During evaluation of In[6]:= {7,28}

During evaluation of In[6]:= {8,41}

During evaluation of In[6]:= {9,59}

During evaluation of In[6]:= {10,83}

During evaluation of In[6]:= {11,117}

During evaluation of In[6]:= {12,161}

During evaluation of In[6]:= {13,218}

During evaluation of In[6]:= {14,294}

During evaluation of In[6]:= {15,397}

During evaluation of In[6]:= {16,530}

During evaluation of In[6]:= {17,704}

During evaluation of In[6]:= {18,936}

During evaluation of In[6]:= {19,1241}

During evaluation of In[6]:= {20,1639}

During evaluation of In[6]:= {21,2169}

During evaluation of In[6]:= {22,2860}

During evaluation of In[6]:= {23,3763}

During evaluation of In[6]:= {24,4935}

During evaluation of In[6]:= {25,6468}

During evaluation of In[6]:= {26,8450}

During evaluation of In[6]:= {27,11031}

During evaluation of In[6]:= {28,14401}

During evaluation of In[6]:= {29,18805}

During evaluation of In[6]:= {30,24542}

During evaluation of In[6]:= {31,32019}

During evaluation of In[6]:= {32,41760}

During evaluation of In[6]:= {33,54447}

During evaluation of In[6]:= {34,70993}

During evaluation of In[6]:= {35,92579}

During evaluation of In[6]:= {36,120670}

During evaluation of In[6]:= {37,157256}

During evaluation of In[6]:= {38,204881}

During evaluation of In[6]:= {39,266915}

During evaluation of In[6]:= {40,347656}

During evaluation of In[6]:= {41,452767}

During evaluation of In[6]:= {42,589626}

During evaluation of In[6]:= {43,767878}

During evaluation of In[6]:= {44,999923}

During evaluation of In[6]:= {45,1302069}

During evaluation of In[6]:= {46,1695354}

During evaluation of In[6]:= {47,2207435}

During evaluation of In[6]:= {48,2874102}

During evaluation of In[6]:= {49,3742053}

During evaluation of In[6]:= {50,4871874}

During evaluation of In[6]:= {51,6342700}

During evaluation of In[6]:= {52,8257432}

During evaluation of In[6]:= {53,10750196}

During evaluation of In[6]:= {54,13995260}

During evaluation of In[6]:= {55,18219685}

Out[6]= {170.046529, Null} *)


--- end edit #3 ---

--- edit #4 ---

I should mention some issues which I don't have time right now to address in the code above.

(1) My square free test is far from optimal. I should have compared all neighbors, then all neighbor pairs, then triples, etc. That would be O(n^2) with n the word length. What I did, in contrast, can repeat character pair comparisons. It might be O(n^2) on average, but it can be cubic and certainly is not the best way to go even for average cases.

(2) The enumeration should use only the "tail" test, that is, only check neighbor sets that contain the new element. This increases speed substantially. Each child test is still O(n^2) worst case though.

(3) There may be ways to improve on (2). For example, if you know that the parent word never required a test of length k neighbors to go all the way to the last elements in each pair member, for all k

(4) Memory remains a serious issue. One way to address that is to pack words using FromDigits[...,3]. This would give a compression ratio of around 20. If one wants to use Compile throughout it gets slightly trickier since we'd want to use an array of machine ints rather than a bigint. If the above spped improvements are useful then maybe one need not bother so much with Compile.

--- end edit #4 ---

-
Correct, I want to count them. –  M.R. Aug 24 '12 at 9:56
Enumerating can be done roughly the same way. At each step where you add an element to an existing sq-free word, you have to check all subwords up to floor(half length of list) that end at this new element with the corresponding subwords to the left. –  Daniel Lichtblau Aug 24 '12 at 14:28
Perhaps if you add , RuntimeOptions -> "Speed", CompilationTarget -> "C" to your notSqrFreeAtTail :D –  belisarius is forth Aug 27 '12 at 3:42
@Verde Thanks, good idea. I thought that would not help because "InlineExternalDefinitions" should just slurp in the definition (or so I thought). Also adding those options gave an ominous error message "Compile::nogen: A library could not be generated from the compiled function. >>" But it works anyway, and around 5x faster than what I had. Fascinating. (And a nice addition to my level of ignorance). –  Daniel Lichtblau Aug 27 '12 at 14:19
The difficulty in enumeration is memory, the list becomes so long that Mathematica dies around n=63. –  M.R. Aug 27 '12 at 19:05

A very simple and straightforward test for square-freeness (and should be reasonably fast) is:

squareFreeQ[str_] := StringFreeQ[str, x__ ~~ x__]


squareFreeQ["0101"]
(* False*)
squareFreeQ["0102012021"]
(* True *)


You can then possibly restrict this further to operate only on certain alphabets using Repeated and Alternatives.

Here's the timings I get for the test string provided in the question (TimeConstrained is used because the code in the question takes forever... still running after 600 seconds!):

TimeConstrained[squareFreeQ@benchmark // Timing, 20, "Out of time"]
(* {1.02648, True} *)

TimeConstrained[SquareFree@benchmark // Timing, 20, "Out of time"]
(* "Out of time!" *)

-
+1, very nice use of string patterns. –  Leonid Shifrin Aug 20 '12 at 20:39
But is it faster than my code? I'm running it now. –  M.R. Aug 20 '12 at 20:41
@Mike Yes, it should be. Mine ran in about 1 second for your test data, whereas yours is still running... –  R. M. Aug 20 '12 at 20:42
LOL, I should have started yours first! +1 –  M.R. Aug 20 '12 at 20:43
@Mike String patterns in Mathematica usually get compiled down to regexps, AFAIK, in a rather wide class of cases, so they are different from normal Mathematica patterns and can be much more efficient (even though using similar syntax). –  Leonid Shifrin Aug 20 '12 at 22:15

If you need the ultimate speed, the following compiled code will be about 20 - 30 times faster than the elegant string-pattern based solution of @R.M. (but, of course, as many times longer and uglier):

With[{part = CompileGetElement},
squareFreeQLSC =
Compile[{{ll, _Integer, 1}},
Module[{res = 0, ctr = 1, sctr = 1, len = 0, start = 0, i = 0},
For[i = 0, i < Length[ll], i++,
If[res != 0, Break[]];
len = Length[ll] - i;
start = i;
ctr = 2;
If[len > 1 && part[ll, start + 1] == part[ll, start + 2],
res = 1;
Break[]
];
While[ctr < len,
If[res != 0, Break[]];
sctr = 1;
While[ ctr <= len && part[ll, ctr + start] != part[ll, start + 1],
ctr++
];
If[ctr > len/2 + 1, Break[]];
While[ctr < len  &&
part[ll, ++ctr + start] == part[ll, sctr + start + 1],
sctr++;
If[ctr == 2*sctr,
res = 1;
Break[]
];
]
]
]; (* For *)
res
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]; (* With *)


The above code defines a function which detects square lists of integers. To adopt it to our needs, we can use the following function:

ClearAll[squareFreeQLS];
squareFreeQLS[str_String] :=
squareFreeQLSC[ToCharacterCode[str]] == 0;


We use it as:

squareFreeQLS[benchmark]

(* True *)

squareFreeQLS[benchmark <> "0120121"]

(* False *)


Here are the timings on my machine:

Do[squareFreeQ[benchmark], {10}] // AbsoluteTiming
Do[squareFreeQLS[benchmark], {10}] // AbsoluteTiming

(*
{15.1826172, Null}
{0.5966797, Null}
*)


Despite having this 20-30 times speedup, I must confess that I would probably go with @R.M.'s solution, given its elegance, simplicity, and also development time (it took me a while to debug and optimize this one), unless you really need the ultimate speed. Of course, there may be more efficient algorithms than the rather straight-forward one I used here, which may be another source of speed-up.

-
Brilliant code, thanks Leonid! –  M.R. Aug 21 '12 at 3:38
@Mike My pleasure :). Looks like I drifted towards Compile in my recent answers. –  Leonid Shifrin Aug 21 '12 at 3:53

Edit

This one is much simpler than those I posted before . And very efficient

Timing@StringFreeQ[benchmark, RegularExpression["(.+)\\1"]]


Previous posts:

Timings done on a VERY slow machine:

Timing@Not@StringMatchQ[benchmark, RegularExpression[".*(.{1,1000})\\1(.*)"]]
(*
-> {0.735, True}
*)


Edit

There is a problem if the repeated string has more than 1325 chars long, because a Max recursion error arises. \$RecursionLimit does not work. So, we could do (still reasonably fast:

mod = 1325;
Timing[Nor @@
Table[StringMatchQ[benchmark, RegularExpression[
".*(.{" <> ToString[1 + mod i] <> "," <> ToString[mod (i + 1)] <>"})\\1(.*)"]],
{i, 0, IntegerPart[StringLength@benchmark/2/mod]}]]
(*
--- > {3., True}
*)


@R.M's string pattern function is 25% slower here. Anyone can confirm?

Edit Breaking down the regexp

".*(.{1,1000})\\1(.*)"

.*             --> Any char (.) zero or more times (*)
(.{1,1000})    --> Any char (not necessarily the same) from 1 to 1000 times
\\1            --> The result of the first parenthesized expression again ...
(.*)           --> Any char (.) zero or more times (*)


The \\1 is the token that specifies the repetition.

-
Yes, I can confirm it is true. Could you also break down the regex and explain? Btw, you need to add a Not to get the OP's definition –  R. M. Aug 20 '12 at 22:36
@R.M I explained the first one. The second is the same with a variable {1,1000} thingy –  belisarius is forth Aug 20 '12 at 22:44
I used squareFreeQ[str_] := StringFreeQ[str, RegularExpression["(\\d+)\\1"]] myself... –  J. M. is back. Aug 21 '12 at 0:34
@J.M. I was going to post Timing@StringFreeQ[benchmark, RegularExpression["(.+)\\1"]]` –  belisarius is forth Aug 21 '12 at 0:39