There are many concrete examples about how to solve a specific system of linear congruences with two variables, for example Problem with solving systems of linear congruence of two variables and How do you solve linear congruences with two variables.. I'd like to ask how to solve a general system of linear congruences with two variables, because it seems difficult to find practical algorithms to solve it :(
$\left\{ \begin{eqnarray} a_1x+b_1y&=c_1 \ (mod \ m)\\ a_2x+b_2y&=c_2 \ (mod \ m)\\ \cdots&\\ a_nx+b_ny&=c_n \ (mod \ m)\\ \end{eqnarray} \right. $
Update: For the case $n=2$, I use Smith normal form to rewrite the above into matrix form $ \begin{pmatrix} a_1 & b_1\\ a_2 & b_2\\ \end{pmatrix}=S*\begin{pmatrix} gcd(a_1,b_1,a_2,b_2) & 0\\ 0 & \frac{|a_1b_2-a_2b_1|}{gcd(a_1,b_1,a_2,b_2)}\\ \end{pmatrix}*T $.
Therefore, $ \begin{pmatrix} a_1 & b_1\\ a_2 & b_2\\ \end{pmatrix}\cdot \begin{pmatrix} x\\ y\\ \end{pmatrix} = \begin{pmatrix} gcd(a_1,b_1,a_2,b_2) & 0\\ 0 & \frac{|a_1b_2-a_2b_1|}{gcd(a_1,b_1,a_2,b_2)}\\ \end{pmatrix}\cdot T\cdot \begin{pmatrix} x\\ y\\ \end{pmatrix} \equiv S^{-1} \begin{pmatrix} c_1\\ c_2\\ \end{pmatrix}=\begin{pmatrix} c_1'\\ c_2'\\ \end{pmatrix} \ (mod \ m) $
Let $a_1'=gcd(a_1,b_1,a_2,b_2)$ and $b_2'=\frac{|a_1b_2-a_2b_1|}{gcd(a_1,b_1,a_2,b_2)}$. Note that $T\cdot \begin{pmatrix} x\\ y\\ \end{pmatrix}=\begin{pmatrix} x'\\ y'\\ \end{pmatrix}$ is one-one mapping w.r.t. $\begin{pmatrix} x\\ y\\ \end{pmatrix}$. Thus the number of solutions is $gcd(a_1,b_1,a_2,b_2,m)\cdot gcd(\frac{|a_1b_2-a_2b_1|}{gcd(a_1,b_1,a_2,b_2)},m)$ if the system is feasible, otherwise $0$. (One corner case is $a_1'=0$ which means one equation holds trivially, in this case if it is feasible, the number of solutions is $gcd(a_1,b_1,a_2,b_2,c)\cdot c$; a general condition is $gcd(a_1',m) \mid c_1'$ and $gcd(b_2',m) \mid c_2'$.)
