How does one solve a system like the following:
$$a_{11}x_{1}+a_{12}x_{2}+...a_{1n}x_{n}=b_{1}\space (mod \space p^{k} )\\\vdots \\ a_{n1}x_{1}+a_{n2}x_{2}+...{a_n}x_{n}=b_{n}\space (mod \space p^{k}) $$ Where $a_{11}, a_{12},\dots ,a_{n(n-1)},a_{nn} $ and $b_{i}$ are $(mod \space p^{k})$ integers?
I could only find linear congruence systems where the Chinese remainder theorem helps but that's not the case here. The only solution that seems good is adding $y_{i} *p^{k}$ to each line and solving the linear system but then we will have $n$ equations for $2n$ uknowns. Is there a better way of solving this?
We consider finding one of solutions of $\mathbf{Ax}=\mathbf{b} \, (\mathrm{mod} p^{k})$ for given $\mathbf{A} \in \mathbb{Z}^{m \times n}$ and $\mathbf{b} \in \mathbb{Z}^{m}$.
Let the Smith normal form (SNF) of $\mathbf{A}$ be $\mathbf{D}=\mathbf{LAR}$, where $L$ and $R$ are some unimodular matrices. Now it is sufficient to solve $\mathbf{Dy} = \mathbf{Lb} \, (\mathrm{mod} p^{k})$ for $\mathbf{y} \in \mathbb{Z}^{m}$. After we obtain $\mathbf{y}$, we easily obtain the solution of the original linear systems as $\mathbf{x} = \mathbf{Ry} \, (\mathrm{mod} p^{k})$.
We write the SNF as $\mathbf{D} = \mathrm{diag}(d_{1}, \cdots, d_{r}, 0, \cdots)$. Then, the $i$th element of $\mathbf{y}$ is determined by solving $$ d_{i} y_{i} = [\mathbf{Lb}]_{i} \quad (\mathrm{mod} \, p^{k}) $$ for $i=1, \cdots, r$. The above can be solved by the extended Euclidean algorithm. If $\mathrm{GCD}(d_{i}, p^{k})$ does not divide $p^{k}$ for some $i$, no solution exists.
Note that multiple solutions $\mathbf{x}$ may exist for $k \geq 2$. If we need all the solutions, the next task is to construct an integer lattice formed by $\mathbf{Ax}=\mathbf{0} \, (\mathrm{mod}\,p^{k})$. I mention this problem can also be solved by a lattice algorithm shown in this lecture note.
The generalization from prime powers $p^{k}$ to other composite numbers $l$ is straightforward. After factoring $l$ and solving the linear equations for each prime power, we can get a solution of $\mathbf{Ax} = \mathbf{b} \, (\mathrm{mod}\, n)$ by the Chinese remainder theorem.
P.S. I have also encountered the situation to solve the modular linear equation and I cannot find a useable implementation to solve it. So, I have implemented the above method in my Python package.