Let $p$ be a prime number and $k$ be a positive integer. How do I determine the number of solutions to a set of equations in variables $0\leq x_1,...,x_n
Let $p$ be a prime number and $k$ be a positive integer. How do I determine the number of solutions to a set of equations in variables $0\leq x_1,...,x_n
You can find not only the number of solutions but also the solutions using Gaussian elimination as follows (working in the ring $\mathbb{Z}/p^k\mathbb{Z}$):
Always select an invertible pivot if possible. The operations in Gaussian elimination don't change the solution set as long as the pivots are invertible. If you encounter a situation in which there's no invertible candidate for the pivot, that means that all candidates contain a factor of $p$. Then you know that the highest digit in the base-$p$ representation of the corresponding variable can be chosen arbitrarily. Then choose a pivot that only has one factor of $p$, and apply elimination like you would "usually", but without multiplying by the common factor of $p$. If there's no pivot with only one factor of $p$, then all candidates have a factor of $p^2$ in common, so you know that the second-highest digit of the corresponding variable can be chosen arbitrarily, and so on.
In this manner, you can perform the entire elimination, always choosing a pivot with the lowest number of factors of $p$, and bring the matrix into row echelon form, which gives you the solution set as usual, keeping in mind that where a pivot isn't invertible, you can choose as many high digits of the corresponding variable arbitrarily as the pivot contains factors of $p$.
In addition to doing invertible row operations which don't change the solution set: that is,
you can also do invertible column operations, which correspond to a change of variable and don't change the number of solutions.
Using invertible row and column operations, any matrix over $\mathbb{Z} / p^k \mathbb{Z}$ can be made into a diagonal matrix following the ideas of computing the Smith normal form, and it is very easy to determine the number of solutions to a diagonal system of equations.
The basic algorithm is as follows:
Note that the Euclidean algorithm is much simpler in $\mathbb{Z} / p^k \mathbb{Z}$ than it is over the integers, because you can divide by anything relatively prime to $p$; e.g. modulo $125$, you can compute $\gcd(10, 75) = \gcd(10/2, 75/3) = \gcd(5,25) = 5$.
Since you only care about the number of solutions, in the third bullet point, after zeroing out the first column the entries in the first row no longer matter, so you could just scribble them out and ignore them for the rest of the calculation.