Complexity of reachability in linear dynamical systems over finite fields Let $A$ be a matrix over the finite field $\mathbb{F}_2 = \{0,1\}$ and $x$, $y$ be vectors of the space $\mathbb{F}_2^n$. I am interested in the computational complexity of deciding whether there exists $t \in \mathbb{N}$ such that $A^t x = y$, i.e., in the reachability problem for linear dynamical systems over finite fields.

The problem is clearly in $\mathbf{NP}$ (guess $0 \le t &lt; 2^n$ and compute $A^t$ in polynomial time by repeated squaring). Me and my colleagues were also able to prove the $\mathbf{NP}$-completeness of the related problem of establishing whether there exists $t \in \mathbb{N}$ such that $A^t x \ge y$, where $\ge$ is componentwise inequality.

This problem seems quite natural, but I have not been able to find references to its computational complexity in the literature, probably because I am not aware of the exact terminology. Do you know if the problem with the equality is $\mathbf{NP}$-complete or if it is actually in $\mathbf{P}$?
For clarity, I'm going to generalize your question to be over characteristic $p
&gt; 0$ (with base field $\mathbb{F}_q$) instead of the specific case of $p=q=2$.
I'll take $p$ and $q$ as a fixed constants;
I'll leave it to the reader to figure out what the exact dependence on these parameters is,
as there are some tradeoffs that can be made.
The end result here is that your problem is roughly equivalent to the discrete
log problem for finite fields of characteristic $p$.

To be more specific,
let the ordinary discrete log problem over extensions of $\mathbb{F}_q$ be,
given an extension field $\mathbb{F}$ of $\mathbb{F}_q$, and $a,b \in \mathbb{F}$,
find any integer $t$ so that $a = b^t$,
or report that none exists.
Let the strong discrete log problem over extensions of $\mathbb{F}_q$ be,
given $\mathbb{F},a,b$ as before, find integers $z,m$ so that
$a = b^t$ for an integer $t$ iff $t = z \pmod{m}$,
or report that no $t$ exists.
Then the following reductions exist:


There is a deterministic mapping reduction from discrete log over extensions of $\mathbb{F}_q$ to your problem.
There is an efficient, deterministic algorithm which solves your problem
when given access to an oracle computing the strong discrete log problem
over extensions of $\mathbb{F}_q$.


Accordingly, I'd consider it unlikely that somebody will post a proof of
$\mathsf{NP}$-hardness or a proof that your problem is in $\mathsf{P}$ in the
near future.

Remark: The strong discrete log problem over extensions of $\mathbb{F}_q$
can be Turing-reduced to the following ostensibly weaker form
(though still seemingly stronger than the ordinary discrete log problem):
Given an extension field $\mathbb{F}$ of $\mathbb{F}_q$, and $a,b \in \mathbb{F}$,
find the least, non-negative integer $t$ so that $a = b^t$.
This follows from the fact that the order of $b$ is one plus the smallest
non-negative $t$ so that $b^{-1} = b^t$.



First reduction:
The claim is that the ordinary discrete log problem over extensions of $\mathbb{F}_{q}$ mapping-reduces to this problem.
This follows the fact that multiplication in $\mathbb{F}_{q^n}$ is a linear transformation
when we view $\mathbb{F}_{q^n}$ as an $n$-dimensional vector space over $\mathbb{F}_q$.
Hence a question of the form $a = b^t$ over $\mathbb{F}_{q^n}$
becomes $\vec{a} = B^t\vec{e}$ over $\mathbb{F}_q$, where $\vec{a},\vec{e}$ are $n$-dimensional vectors,
and $B$ is an $n\times n$ matrix, all over $\mathbb{F}_q$.
The vector $\vec{a}$ can be easily computed from $a$, $B$ from $b$, and
$\vec{e}$ is just the representation of $1 \in \mathbb{F}_{q^n}$, which can be
written down efficiently.
This appears to still be a hard case of the general discrete log problem,
even with $p=q=2$ (but growing $n$, of course).
In particular, people are still competing to see how far out they can compute it.



Second reduction:
The claim is that your problem reduces to the strong discrete log problem over extensions of $\mathbb{F}_q$.
This reduction has a few pieces to it, so forgive the length.
Let the input be the $n$-dimensional vectors $x,y$ and $n\times n$ matrix $A$, all over $\mathbb{F}_q$;
the goal is to find $t$ so that $y = A^tx$.

The basic idea is to write $A$ in Jordan canonical form (JCF),
from which we can reduce testing $y = A^tx$ to the strong discrete log problem
with some straightforward algebra.

One reason for using a canonical form under similarity of matrices
is that if $A = P^{-1}JP$, then $A^t = P^{-1}J^tP$.
Hence we can transform $y = A^tx$ to $(Py) = J^t(Px)$,
where now $J$ is in a much nicer format than the arbitrary $A$.
The JCF is a particularly simple form, which enables the rest of the algorithm.
So from now on, assume that $A$ is already in JCF,
but also allow that $x,y,$ and $A$ may have entries in an extension field of
$\mathbb{F}_q$.

Remark: There are some subtleties that arise from working with the JCF.
Specifically, I'll assume that we can do field operations within any extension
of $\mathbb{F}_q$ (no matter how large) in one time step,
and that we can compute the JCF efficiently.
a priori, this is unrealistic, because working with the JCF may require
working in an extension field (the splitting field of the characteristic polynomial)
of exponential degree.
However, with some care, and using the fact that we're working over a finite field,
we can circumvent these issues.
In particular,
we will associate with each Jordan block a field $\mathbb{F}'$ of degree at most $n$ over $\mathbb{F}_q$
so that all the entries in the Jordan block and the corresponding elements of $x$,$y$
all live within $\mathbb{F}'$.
The field $\mathbb{F}'$ may differ from block to block,
but using this ``mixed representation'' allows for an efficient description of the JCF,
which moreover can be found efficiently.
The algorithm described in the remainder of this section only needs to work
with one block at a time,
so as long as it does its field operations within the associated field $\mathbb{F}'$,
the algorithm will be efficient.
[end remark]

The use of JCF gives us equations of the following form,
with each equation corresponding to a Jordan block:
$$
\begin{bmatrix}
  y_1 \\
  y_2 \\
  y_3 \\
  \vdots \\
  y_{k-1} \\
  y_{k} \\
\end{bmatrix}
=
\begin{bmatrix}
\lambda &amp; 1       &amp;         &amp;        &amp;         &amp;         \\
    &amp; \lambda &amp; 1       &amp;        &amp;         &amp;         \\
    &amp;         &amp; \lambda &amp; 1      &amp;         &amp;         \\
    &amp;         &amp;         &amp; \ddots &amp;         &amp;         \\
    &amp;         &amp;         &amp;        &amp; \lambda &amp; 1       \\
    &amp;         &amp;         &amp;        &amp;         &amp; \lambda \\
\end{bmatrix}^t
\begin{bmatrix}
  x_1 \\
  x_2 \\
  x_3 \\
  \vdots \\
  x_{k-1} \\
  x_{k} \\
\end{bmatrix}$$

The algorithm will handle each block separately.
In the general case, for each block, we'll have a query for our strong discrete log oracle,
from which the oracle will tell us a modularity condition, $t = z \pmod{m}$.
We'll also get a set $S \subseteq \{0,1,\cdots,p-1\}$ so that
$\bigvee_{s\in S}\left[ t = s \pmod{p} \right]$
must hold.
After processing all the blocks,
we'll need to check that there is a choice of $t$ that satisfies the conjunctions of all these conditions.
This can be done by making sure there is a common element $s$ in all the sets $S$
so that the equations $t = s \pmod p$ and $t = z_j \pmod{m_j}$ are all simultaneously satisfied,
where $j$ ranges over the blocks.

There are also some special cases that arise throughout the procedure.
In these cases,
we'll get conditions of the form $t &gt; \ell$ for some value of $\ell$,
or of the form $t = s$ for some specific integer $s$, from certain blocks,
or we might even find that no $t$ can exist.
These can be incorporated into the logic for the general case without issue.

We now describe the subprocedure for handling each Jordan block.
Fix such a block.

Begin by focusing on just the last coordinate in the block.
The condition $y = A^tx$ requires that $y_k = \lambda^t x_k$.
In other words, it's an instance of the discrete log problem in some field extension of $\mathbb{F}_q$.
We then use an oracle to solve it,
which either results in no solution,
or else gives a modularity condition on $t$.
If "no solution" is returned, we return indicating such.
Otherwise,
we get a condition $t = z \pmod{m}$,
which is equivalent to $y_k = \lambda^t x_k$.

To handle the other coordinates,
we start with the following formula
(see, eg, here):
$$
\begin{bmatrix}
\lambda &amp; 1       &amp;         &amp;        &amp;         &amp;         \\
    &amp; \lambda &amp; 1       &amp;        &amp;         &amp;         \\
    &amp;         &amp; \lambda &amp; 1      &amp;         &amp;         \\
    &amp;         &amp;         &amp; \ddots &amp;         &amp;         \\
    &amp;         &amp;         &amp;        &amp; \lambda &amp; 1       \\
    &amp;         &amp;         &amp;        &amp;         &amp; \lambda \\
\end{bmatrix}^t
=
\begin{bmatrix}
\lambda^t &amp; \binom{t}{1}\lambda^{t-1} &amp; \binom{t}{2}\lambda^{t-2} &amp; \cdots &amp; \cdots &amp; \binom{t}{k-1}\lambda^{t-k+1} \\
 &amp; \lambda^t &amp; \binom{t}{1}\lambda^{t-1} &amp; \cdots &amp; \cdots &amp; \binom{t}{k-2}\lambda^{t-k+2} \\
 &amp;  &amp; \ddots &amp; \ddots &amp; \vdots &amp; \vdots\\
 &amp;  &amp; &amp; \ddots &amp; \ddots &amp; \vdots\\
 &amp;  &amp; &amp;  &amp; \lambda^t &amp; \binom{t}{1}\lambda^{t-1}\\
 &amp;  &amp;  &amp;  &amp;  &amp; \lambda^t
\end{bmatrix}$$
First, let's take care of the case in which $x_k = 0$.
Since we already have the modularity condition which implies $y_k = \lambda^t x_k$,
we can assume that $y_k = 0$ also.
But then we can just reduce to focusing on the first $k-1$ entries of $x$ and $y$,
and the top left $(k-1)\times (k-1)$ submatrix of the Jordan block.
So from now on, assume that $x_k \ne 0$.

Second, we'll handle the case in which $\lambda = 0$.
In this case,
the powers of the Jordan block have a special form,
and force either $t = z$ for some $z \le k$,
or else $t &gt; k$, with no other conditions.
I won't belabor the cases, but suffice it to say that each can be checked for efficiently.
(Alternatively, we could reduce to the case where $A$ is invertible; see my comment on the question.)

Finally, we arrive at the general case.
Since we already have the modularity condition which implies that $y_k = \lambda^t x_k$,
we can assume that condition holds, and use $y_k x_k^{-1}$ as a stand-in for $\lambda^t$.
More generally, we can use $y_kx_k^{-1}\lambda^{-z}$ to represent $\lambda^{t-z}$.
Thus we need to check if the following system holds for some choice of $t$:
$$
\begin{bmatrix}
  y_1 \\
  y_2 \\
  y_3 \\
  \vdots \\
  y_{k-1} \\
  y_{k} \\
\end{bmatrix}
=
\begin{bmatrix}
y_kx_k^{-1} &amp; \binom{t}{1}y_kx_k^{-1}\lambda^{-1} &amp; \binom{t}{2}y_kx_k^{-1}\lambda^{-2} &amp; \cdots &amp; \cdots &amp; \binom{t}{k-1}y_kx_k^{-1}\lambda^{-(k-1)} \\
 &amp; y_kx_k^{-1} &amp; \binom{t}{1}y_kx_k^{-1}\lambda^{-1} &amp; \cdots &amp; \cdots &amp; \binom{t}{k-2}y_kx_k^{-1}\lambda^{-(k-2)} \\
 &amp;  &amp; \ddots &amp; \ddots &amp; \vdots &amp; \vdots\\
 &amp;  &amp; &amp; \ddots &amp; \ddots &amp; \vdots\\
 &amp;  &amp; &amp;  &amp; y_kx_k^{-1} &amp; \binom{t}{1}y_kx_k^{-1}\lambda^{-1}\\
 &amp;  &amp;  &amp;  &amp;  &amp; y_kx_k^{-1}
\end{bmatrix}
\begin{bmatrix}
  x_1 \\
  x_2 \\
  x_3 \\
  \vdots \\
  x_{k-1} \\
  x_{k} \\
\end{bmatrix}
$$
Observe that whether the equation holds depends only on $t \pmod{p}$;
this is because the dependence on $t$ is only polynomial,
$t$ must be an integer,
and the above equations are over a field of characteristic $p$.
Hence we can just try each value of $t \in \{0,1,\ldots,p-1\}$ separately.
The set $S$ we will return is just the choices of $t$ for which the system is
satisfied.

So now, except for some special cases, the per-block subprocedure has found a modularity condition $t = a \pmod{m}$,
and a set $S$ so that one of $t = s \pmod{p}$ must hold for some $s \in S$.
These conditions are equivalent to $y = A^tx$ within this specific Jordan block.
So we return these from the subprocedure.
The special cases either conclude that no $t$ can exist (in which case the subprocedure immediately returns an indication of that),
or else we have a modularity condition $t = a \pmod{m}$ and some special condition like $t = s$ for an integer $s$,
or $t &gt; \ell$ for some integer $\ell$.
In any case, the conditions involved are all equivalent to $y = A^tx$ within this Jordan block.
So, as mentioned above, the subprocedure just returns these conditions.

This concludes the specification of the per-block subprocedure,
and of the algorithm as a whole.
It's correctness and efficiency follow from the preceding discussion.



Subtleties with using JCF in second reduction:
As mentioned in the second reduction, there are some subtleties that arise from
working with the JCF.
There are a few observations for mitigating these problems:


Extensions of finite fields are normal.
This means that if $P$ is a an irreducible polynomial over $\mathbb{F}_q$,
then any extension of $\mathbb{F}_q$ containing a root of $P$
contains all the roots of $P$.
In other words, the splitting field of an irreducible polynomial $P$
of degree $d$ has degree only $d$ over $\mathbb{F}_q$.
There is a generalization of the Jordan canonical form,
called the primary rational canonical form (PRCF),
which does not require field extensions to be written down.
In particular, if $A$ is a matrix with entries in $\mathbb{F}_q$,
then we can write $A = P^{-1}QP$ for some matrices $P,Q$ with entries in $\mathbb{F}_q$,
where moreover $Q$ is in PRCF.
Additionally, if we pretend that the entries $A$ live in a field $\mathbb{F}'$
extending $\mathbb{F}_q$ which contains all the eigenvalues of $A$,
then $Q$ will in fact be in JCF.
Thus we can view computing the JCF of $A$ as a special case of computing the PRCF.
Using the form of the PRCF,
we can factor computing the JCF of $A$ as


computing the PRCF of $A$ over $\mathbb{F}_q$
computing the PRCF of each block $C$ (borrowing the notation from the Wikipedia article) in the PRCF of $A$,
over an extension field $\mathbb{F}'$,
where $\mathbb{F}'$ is chosen to contain all the eigenvalues of $C$


The key advantage with this factorization is that the characteristic
polynomials of the blocks $C$ will all be irreducible,
and hence, by our first observation,
we can choose $\mathbb{F}'$ to have degree the size of $C$ (which is at most $n$) over $\mathbb{F}_q$.
The downside is that now we have to use different extension fields to
represent each block of the JCF, so the representation is atypical and complicated.


Thus, given the ability to compute the PRCF efficiently,
we can compute a suitable encoding of the JCF efficiently,
and this encoding is so that working with any particular block of the JCF can be done within
an extension field of degree at most $n$ over $\mathbb{F}_n$.

As for computing the PRCF efficiently,
the paper
"A Rational Canonical Form Algorithm" (K. R. Matthews, Math. Bohemica 117 (1992) 315-324)
gives an efficient algorithm to compute the PRCF when the factorization of the characteristic
polynomial of $A$ is known.
For fixed characteristic (such as we have),
factoring univariate polynomials over finite fields can be done in deterministic polynomial time
(see eg "On a New Factorization Algorithm for Polynomials over Finite Fields" (H. Niederreitter and R. Gottfert, Math. of Computation 64 (1995) 347-353).),
so the PRCF can be computed efficiently.