
MATHEMATICA

g = Function[{sq, p}, Module[{l = Length[sq]},
Do[If[sq[[i]] == sq[[j]], Return[p^(j  1)  p^(i  1)]],
{j, 2, l}, {i, 1, j  1}]]];
MPM = Algebra`MatrixPowerMod;
EventualPeriod = Function[{m, v, p},
Module[{n = Length[m], w, sq, k, primes},
sq = NestList[(MPM[#, p, p]) &, m, n];
w = Mod[Last[sq].v, p];
sq = Map[(Mod[#.w, p]) &, sq];
k = g[sq, p];
If[k == Null, k = p^n Apply[LCM, Table[p^r  1, {r, 1, n}]]];
primes = Map[First, FactorInteger[k]];
primes = Select[primes, (# > 1) &];
While[Length[primes] > 0,
primes = Select[primes, (Mod[k, #] == 0) &];
primes = Select[primes, (Mod[MPM[m, k/#, p].w, p] == w) &];
k = k/Fold[Times, 1, primes];
]; k ]];
mat = Function[{n}, Table[Boole[Abs[i  j] == 1], {i, 1, n}, {j, 1, n}]];
vec = Function[{n}, Table[Boole[i == 1], {i, 1, n}]];
Table[EventualPeriod[mat[2 n], vec[2 n], 2], {n, 1, 100}]
