  We show that the density $\mu$ of the Smith normal form (SNF) of a random
integer matrix exists and equals a product of densities $\mu_{p^s}$ of SNF over
$\mathbb{Z}/p^s\mathbb{Z}$ with $p$ a prime and $s$ some positive integer. Our
approach is to connect the SNF of a matrix with the greatest common divisors
(gcds) of certain polynomials of matrix entries, and develop the theory of
multi-gcd distribution of polynomial values at a random integer vector. We also
derive a formula for $\mu_{p^s}$ and compute the density $\mu$ for several
interesting types of sets. Finally, we determine the maximum and minimum of
$\mu_{p^s}$ and establish its monotonicity properties and limiting behaviors.
