Cycle detection with Mathematica

Question

Suppose we have cellular automaton on a network. For simplicity, we will use matrix notation.

ClearAll[adjMatrix, initStates, nodeStep, allStep];
(*Adjacency matrix*)
adjMatrix = {{0, 1, 1, 0}, {0, 1, 0, 0}, {1, 0, 1, 1}, {0, 0, 1, 0}};
(*Initial states of nodes*)
initStates ={0, 1, 1, 0};  

nodeStep[adjMatrix_, states_, node_] :=
  With[{inputs = Pick[states, adjMatrix[[All, node]], 1]},
    (*Any suitable function here*)
    BitXor @@ inputs
   ];

allStep[adjMatrix_, states_] := 
  nodeStep[adjMatrix, states, #] & /@ Range[Length@states];

Starting from some initial state, the function allStep is applied iteratively. It is known that sooner or later we will get a cycle.

For test example:
{0, 1, 1, 0} → {1, 1, 1, 1} → {1, 0, 1, 1} → {1, 1, 1, 1} → ...
(period 2)

But first, not necessarily straight from the initial state. Second, about the cycle length (period), it is only known that it is smaller than $2^{size}$

I have not been able to find a way to simultaneously detect the cycle and determine its length. For detection I use

data = NestWhileList[allStep[adjMatrix, #]&, initStates, Unequal, All];

and then we can find the length of the cycle.
Several ways have been suggested here.

• FindRepeat. It fails with period 1
• SequencePosition. I do not understand what should be M in my case: SequencePosition[data, Take[data, M]]?
• FindTransientRepeat. It works, but much slower than brutal force method
• First@Differences@Flatten@Position[data, Last@data]

Some timings:

data = ContinuedFraction[(Sqrt[12] + 2)/7, 100004];
Timing[Length@Last@FindTransientRepeat[data, 2]]
{0.499203, 6}
Timing[r = SequencePosition[data, Take[data, -10]]; 
 r[[-1, 1]] - r[[-2, 1]]]
{0.0156001, 6}
Timing[First@Differences@Flatten@Position[data, Last@data]]
{0.0468003, 6}

Answer 1

The best solution is Brent's algorithm.

brentCycleDetection[adjMatrix_, states_] :=
  Module[{power = 1, lam = 1, tortoise = states, 
    hare = allStep[adjMatrix, states]},
   While[tortoise != hare,
     If[ power == lam,
                 tortoise = hare;
                 power *= 2;
                 lam = 0];
    hare = allStep[adjMatrix, hare];
    lam += 1;
    ];
   lam
   ];

Here is typical example:

size = 13; adjMatrix = RandomInteger[1, {size, size}]; states = 
 RandomInteger[1, size];

Timing[brentCycleDetection[adjMatrix, states]]
{0.140401, 510}

Timing[myCycleDetection[adjMatrix, states]]
{1.54441, 510}