We show that the density μ of the Smith normal form (SNF) of a random integer matrix exists and equals a product of densities μps of SNF over Z/psZ 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 μps and determine the density μ for several interesting types of sets.