# 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. Aug 22, 2012 at 22:09
– M.R.
Aug 22, 2012 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, 2012 at 22:39
• @Rojo, this shows that the number will be something around 10^13. Aug 23, 2012 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, 2012 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, 2012 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. Aug 24, 2012 at 14:28
• Perhaps if you add , RuntimeOptions -> "Speed", CompilationTarget -> "C" to your notSqrFreeAtTail :D Aug 27, 2012 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). Aug 27, 2012 at 14:19
• The difficulty in enumeration is memory, the list becomes so long that Mathematica dies around n=63.
– M.R.
Aug 27, 2012 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. Aug 20, 2012 at 20:39
• But is it faster than my code? I'm running it now.
– M.R.
Aug 20, 2012 at 20:41
• @Mike Yes, it should be. Mine ran in about 1 second for your test data, whereas yours is still running...
– rm -rf
Aug 20, 2012 at 20:42
• LOL, I should have started yours first! +1
– M.R.
Aug 20, 2012 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). Aug 20, 2012 at 22:15

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
– rm -rf
Aug 20, 2012 at 22:36
• @R.M I explained the first one. The second is the same with a variable {1,1000} thingy Aug 20, 2012 at 22:44
• I used squareFreeQ[str_] := StringFreeQ[str, RegularExpression["(\\d+)\\1"]] myself... Aug 21, 2012 at 0:34
• @J.M. I was going to post Timing@StringFreeQ[benchmark, RegularExpression["(.+)\\1"]] Aug 21, 2012 at 0:39

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, 2012 at 3:38
• @Mike My pleasure :). Looks like I drifted towards Compile` in my recent answers. Aug 21, 2012 at 3:53