I have a system of linear equations over a finite field $\mathbb F_p \cong \mathbb Z_p$, and I'm interested in the decision problem of whether there exists a solution where all of the variables $x_i$ are in the set $\{0, 1\} \subset \mathbb F$. In particular, I'm trying to determine whether this problem is NP-hard.
Example
One system of equations over $\mathbb F_3$ is: $$ \begin{alignat*}{2} &x_1\begin{bmatrix}1 \\ 1 \\ 0\end{bmatrix} + x_2\begin{bmatrix}0 \\ 0 \\ 0\end{bmatrix} + x_3\begin{bmatrix}0 \\ 2 \\ 0\end{bmatrix} + x_4\begin{bmatrix}1 \\ 0 \\ 1\end{bmatrix} \\ + \, &x_5\begin{bmatrix}2 \\ 1 \\ 0\end{bmatrix} + x_6\begin{bmatrix}0 \\ 1 \\ 1\end{bmatrix} + x_7\begin{bmatrix}1 \\ 0 \\ 0\end{bmatrix} + x_8\begin{bmatrix}1 \\ 2 \\ 0\end{bmatrix} = \begin{bmatrix}0 \\ 0 \\ 2\end{bmatrix}. \end{alignat*} $$ This system of equations is satisfiable with entries in $\{0,1\}^8 \subset \mathbb F^8$, namely $$ \begin{align*} (x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8) &= (0,0,0,1,0,1,1,1) \hspace{1em}\text{or}\\ (x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8) &= (0,1,0,1,0,1,1,1). \end{align*} $$
An unhelpful (?) reduction
One suggestion that was given to me was turning this into a system of quadratic equations in the following way: define auxiliary functions coordinatewise $$ \begin{align*} w_1 &= x_1 + x_4 + 2x_5 + x_7 + x_8 \\ w_2 &= x_1 + 2x_3 + x_5 + x_6 + 2x_8 \\ w_3 &= x_4 + x_6 - 2, \end{align*} $$ and use these to solve the system of quadratic and linear equations $$ w_1 = w_2 = w_3 = x_1^2 - x_1 = \cdots = x_8^2 - x_8 = 0. $$
However, the MQ-problem (Multivariate Quadratic equations over a finite field) is NP-hard, so this reduction doesn't help. However, this set-up is a quite special case, so I'm holding out some hope that the original problem might still be in P.
Is there a polynomial-time algorithm for determining the existence of a solution of linear equations over a finite field with restricted variables? Or is it known if this problem is NP-hard like the MQ-problem?
