### Solving a system of linear congruences using the Smith normal form (SNF) of an integer matrix

A is a nonzero m × n matrix of integers, B is m × 1, X is n × 1. We solve the system of linear congruences AX ≡ B (mod q):

a11x1 + ⋯ + a1nxn ≡ b1 (mod q)
.
.
.
am1x1 + ⋯ + amnxn ≡ bm (mod q)

using the Smith normal form of the integer matrix A.

P and Q are unimodular matrices such that PAQ = diag(d1,...,dr,0,...,0), where r = rank(A) and d1,...,dr are positive integers such that di divides di+1 for 1 ≤ i ≤ r-1.
Write X = QY and K = PB. Then AX ≡ B (mod q) is equivalent to the system of congrences

d1y1 ≡ k1 (mod q)
.
.
.
dryr ≡ kr (mod q)
0 ≡ kr+1 (mod q)
.
.
.
0 ≡ kn (mod q).

Assuming that the first r congruences are soluble (mod q) with solution arrays Y1,...,Yr of lengths c1,...,cr, where ci = gcd(di,q) and that the last n-r congruences are satisfied, we get the complete solution of the form
X = QY, where Y = [y1,...,yr,yr+1,...,yn]t, where yi = Yi,zi and yr+1,...,yn are arbitrary (mod q).

If r = n, we get c1···cr solutions (mod q), otherwise we get c1···crqn-r solutions (mod q).

The augmented matrix [A|B] can be entered either
(i) as a string of m(n+1) integers separated by spaces, or
(ii) cut and pasted from a text file, with entries separated by spaces and each row ended by a newline.

We assume that [A|B] is not the zero matrix (mod q).

The author is grateful to Alan Offer for programming assistance with the recursive construction of the cartesian product.

Enter m, the number of equations (≤ 50):
Enter n, the number of unknowns (≤ 50):
Enter the modulus q (> 1):
Enter the m × (n+1) augmented matrix :