- Let $A\in M_{m,n}(\mathbb{Z})$; we consider an equation in the form $AX=0$ where $X\in \mathbb{Z}^n$. Let $r=rank(A)$.
We write the Smith Form of $A$: $A=P^{-1}SQ^{-1}$ where $S$ is pseudo-diagonal and $P,Q\in GL(\mathbb{Z})$. Then $SQ^{-1}X=0$ or $SY=0$ with $X=QY$. We obtain, up to order, $Y=[0_r,y_1,\cdots,y_{n-r}]^T$ where the $(y_j)$ are the parameters of the general solution. Finally, the $(x_i)_{i\leq n}$ are known linear combinations of the $(y_j)$ and the complexity of the calculation is polynomial.
Each condition $x_k\not= x_l$ is equivalent to a forbidden linear relation $R_{k,l}=0$ between the $(y_j)$. Using -for example- the "isolve" command of Maple, we can obtain a parameterization of the solutions of $R_{k,l}=0$ in the form $y_j=f_j(z_1,\cdots,z_{n-r-1})$. Yet, such a parameterization does not seem very useful because the parameters change with the relations $R_{k,l}$.
- $\mathbb{Z}/N\mathbb{Z}$ is not a PID, except when $N$ is a prime. We consider the case when $N=p_1\cdots p_k$ where the $(p_i)$ are distinct primes. Then it suffices to solve the system over each $\mathbb{Z}/p_i\mathbb{Z}$. In the sequel, we assume that $N=p$ is a prime.
Of course, the set of solutions is finite. I don't know if you want the list of all solutions or just the number of solutions... for example
let $p=17,n=10,r=4$. A test gives the following results: without any inequations, there are $17^6=24,137,569$ solutions. With $2$ inequations, there are $21,381,376$ solutions.
Since $K$ is a field, it is easy to calculate the solutions without any inequations. We may choose the $n-r$ parameters among the $(x_i)$ and the other $(x_i)$ are linear combinations of the previous ones.
About the inequations $x_k\not= x_l$, if we don't want to be satisfied with forbidden equations, then we can transform them into equations in 2 different ways:
$\bullet$ $(x_k-x_l)^{p-1}=1$. We add an equation of degree $p-1$ and we don't add any unknown.
$\bullet$. $(x_k-x_l)k=1$. We add an equation of degree $2$ and an unknown $k$.
We can solve the system using the Grobner basis theory; however, the complexity can become very large with $n, p$.
EDIT. Answer to joro
No, the case when $p=3$ works correctly. For example, let $n=20,r=15$.
We use Grobner's theory and the method $(x_k-x_l)^{p-1}=1$; a test with $6$ inequations gives the following result
The list of $32$ solutions is displayed in $8$ seconds.
If we only want to know the number of solutions, then the result is displayed in $0.2$ second.
If $p=5$ or more generally if $p>3$, then the calculation is longer but is feasible if the number of solutions is of the order of a few tens.
.