Find the order of a class of finite matrices over finite fields Consider matrices $M$ of size $L\times L$ over a finite field $\mathbb{Z}_p$, for simplicity focus on $p$ prime. The size $L$ is even. We want to find the order of a specific class of matrices, namely we want to find the smallest non-zero integer $n$ such that $$M^n=1,$$ where 1 is here the identity matrix of size $L\times L$. The matrix $M$ has a specific block structure, which originates from a specific linear cellular automaton. We have $$M=AB,$$ where
$$A=
  \begin{pmatrix}
    1 & 1 & & && & \\
    1 & -1 && & & & \\
   & &    1 & 1 && &\\
   & &  1 & -1 && & \\
   &&&& \ddots & & \\
   &  & & &&1 & 1\\
     &  & & &&1 & -1\\
  \end{pmatrix}$$
and
$$
  B=
   \begin{pmatrix}
 -1 && & & && 1&\\
   &    1 & 1 && && \\
   &  1 & -1 && && \\
  &&& \ddots & & & \\
    & & &&1 & 1&\\
    & & &&1 & -1&\\
    1 && & & && 1\\
  \end{pmatrix}$$
You can see $B$ as $CAC^{-1}$ where $C$ is a cyclic shift over $L$ variables. Separately $A$ and $B$ have simple properties, but their alternating product becomes complicated.
I did some numerical testing, choosing first $\mathbb{Z}_3$. The interesting thing is (which also motivates me to look deeper into this) that the order $n$ seems to have a complicated behaviour as a function of $L$ and it is not clear to me what is the precise source of this. Sometimes $n$ is very big, seemingly exponentially growing with $L$. For example $n(L=46)=354292=2^2\cdot 23\cdot 3851$. Or $n(L=58)=9565940=2^2\cdot 5\cdot 29\cdot 16493$. But if $L$ is divisible by 6 then we get much lower numbers, $n(60)=120$ for example. I understand that it is natural  to see the prime $p$ somehow reflected in the function $n(L)$, but what I don't understand is how the big primes mentioned above can enter the game.
Update: Further numerical experiments showed that the smallest orbit happens when $L=2p^m$, in which case $n=2L$, we have simple linear growth. It gets complicated and quickly growing for most other $L$.
The origin of the problem is the following cellular automaton. Consider $L$ variables $s_j\in \mathbb{Z}_p$, with  $j=1,2,\dots,L$. We define a dynamics on the model, such that we perform two operations cyclically. In each operation we group two neighbouring variables into a pair, and we alternate the ways how we make the pairing. For each pairing we perform the update
$$(s_j,s_{j+1})\to (s_j+s_{j+1},s_j-s_{j+1})$$
What changes is that at every second update $j$ is even or odd. If you represent this linear operation with matrices, you get $A$ and $B$ and we want to iterate the product $M=AB$.
 A: Big primes enter the picture as follows: Compute the characteristic polynomial $P_L$  of your square matrix $M$ of size $L$.
Factor $P_L$ over your working prime $p$. This gives you the eigenvalues of the involved Jordan-blocs over $\mathbb F_p$. These eigenvalues are elements of field extensions of degree at most $L$ of your groundfield $\mathbb F_p$. You can thus get large primes dividing the order of $M$ if they are divisors of $p^k-1$ for $k$
the degree of (at least) one irreducible polynomial
(over $\mathbb F_p$) dividing $P_L$.
A fast way for computing the order of $M$ is thus to compute the characteristic polynomial $P_L$ of $M$,
factor it over $\mathbb F_p$ and check then if prime-divisors of $p^k-1$ (for $k$ the degree of an involved irreducible polynomial) divide the order.
The order of $M$ over $\mathbb F_p$ is necessarily a divisor of $U=p^{L-1}\prod_{k=1}^{L}(p^k-1)$. One can of course remove a factor $p^k-1$ from this product if $P_L$ mod $p$ has no irreducible divisor modulo $p$ of degree $k$. One can further reduce the size of $U$ by taking greatest common divisors of all involved factors. The exponent $L-1$ of $p$ in $U$ can of course be replaced by $\mu-1$ where $\mu$ is largest occuring multiplicity among the irreducible divisors of $P_L \pmod p$. (There seem indeed to be often multiplicities in $P_L$: This leads to possibly non-trivial Jordan-blocks in $M$ which contribute a factor $p^d$ to the order of $M$ where $d+1$ is the dimension of the largest involved Jordan bloc.)
The computational bottleneck of this approach is the factorization of $p^k-1$. The rest is computationally easy using fast exponentiation. Concretetely,
You have to check if $M^{U/q}$ is still the identity for $q$ any prime divisor of $U$. If this is the case, replace $U$ by $U/q$ (and recheck if you can remove one more power of $q$ from $U$). Going through all prime divisors of $U$ gives you at the end the exact order of $M$ over $\mathbb F_p$.
Experimental observation The characteristic polynomial $P_L$ of the square
matrix $M=AB$ of even size $L$ seems to be given by
$$P_L=(x^2+4)Q^2_{L/2}$$
where $Q_1=1,\quad Q_2=x+2$ and
$Q_n=(x+2)Q_{n-1}-xQ_{n-2}$ for $n\geq 3$.
(The polynomials $Q_n$ are 'almost' orthogonal polynomials and seem to satisfy
the divisibility property $P_a\vert P_b$ if $a\vert b$.) This formula holds for $L\leq 60$ even.
A: Denote $L=2\ell$, and let me assume that $p$ is an odd prime, not dividing $\ell$. Then there are $\ell$ distinct $\ell$-th roots of unity in a suitable extension of ${\mathbb F}_p$.
The Question can be solved by Fourier analysis over ${\mathbb F}_p^\ell$ (which turns out to be a discrete Fourier analysis). Let me split an $L$-vector $x$ into its odd and even parts:
$$y_i=x_{2i-1},\qquad z_i=x_{2i}.$$
Then $x':=Ax$ and $x'':=Bx$ are given by
$$y'=y+z,\quad z'=y-z,\qquad y_i''=-y_i+z_{i-1},\quad z_i''=z_i+y_{i+1}.$$
Let me define the Fourier transform of an $\ell$-vector $v$ by
$$\hat v(\omega)=\sum_iv_i\omega^{-i},$$
where the argument $\omega$ takes values in the group of $\ell$-th roots of unity. The formula above yield
$$\hat y'=\hat y+\hat z,\quad\hat z'=\hat y-\hat z,\qquad\hat y''(\omega)=-\hat y(\omega)+\omega^{-1}\hat z(\omega),\quad\hat z''(\omega)=\hat z(\omega)+\omega\hat y(\omega).$$
We infer that $X:=Mx$ is given by
$$X(\omega)=\hat M(\omega)x(\omega),\qquad\hat M(\omega):=\begin{pmatrix} \omega-1 & 1+\omega^{-1} \\ -1-\omega & \omega^{-1}-1 \end{pmatrix}.$$
The characteristic polynomial of $\hat M(\omega)$ is
$$P_\omega(\lambda)=\lambda^2-(\omega+\omega^{-1}-2)\lambda+4.$$
The spectrum of $M$ is thus the union of the pairs $(\lambda_\omega,\mu_\omega)$ of roots of the $P_\omega$'s.
The blocks $M(\omega)$ are diagonalisable, unless an $\ell$-root of unity satisfies $\omega^2-6\omega+1=0$.
The order of $M$ in ${\bf M}_L({\mathbb F}_p)$ is the lcm of the orders of the eigenvalues $\lambda_\omega,\mu_\omega$ in the algebraic closure (actually some finite extension),  mulitiplied by $2$ in the exceptional case that this lcm is odd (unlikely) and the order of the solutions of $t^2-6t+1=0$ divide $\ell$.
In other words, defining the vectors $v_\omega=(1,\omega,\ldots,\omega^{\ell-1})$, we see that the subspaces defined by $y,z\in{\rm Span}(v_\omega)$ are stable under both $A$ and $B$, and we are able to compute the spectrum of $AB$ by studying its restriction to these $2$-dimensional spaces. Notice that the action on the spaces associated with $v_{\omega}$ and $v_{\omega^{-1}}$ are conjugate to each other.
