The computation requires two processes:
- row operations of the type used in Gaussian elimination (with some restrictions because we require that the integer equivalence class be preserved) and the corresponding column operations,
- the Euclidean algorithm for finding the greatest common divisor of two integers.
Specifically, you are allowed to
- interchange two rows or two columns,
- multiply a row or column by $\pm1$ (which are the invertible elements in $\mathbf{Z}$),
- add an integer multiple of row to another row (or an integer multiple of a column to another column).
The first goal is to reach diagonal form. Let's first work on column 1: using operation 3 for rows repeatedly, you can, by following the Euclidean algorithm, form a row whose first element is the GCD of the elements in column 1. You can then obtain a matrix with the GCD in the $(1,1)$ position and zeroes in the rest of column 1.
Now work on row 1: do the same thing, but using column operations; eventually you will have the GCD of row 1 in the $(1,1)$ position, and zeroes elsewhere in row 1. You will most likely have messed up column 1, but that's OK. Go back and redo column 1, then redo row 1, and repeat until all elements in row and column 1 are 0 except for the $(1,1)$ element. This process is guaranteed to terminate because the GCD gets smaller each time.
Now you can move on to row/column 2, and repeat the process. Continue for row/column 3, and so on, until you have reached diagonal form.
You may not be done at this point, because the diagonal elements may not satisfy the divisibility requirement of the Smith normal form. You can, however, enforce this by some additional moves as in the following example:
$$
\begin{bmatrix}8 & 0\\0 & 12\end{bmatrix}\longrightarrow\begin{bmatrix}8 & 8\\0 & 12\end{bmatrix}\longrightarrow\begin{bmatrix}8 & 8\\-8 & 4\end{bmatrix}\longrightarrow\begin{bmatrix}24 & 8\\0 & 4\end{bmatrix}\longrightarrow\begin{bmatrix}24 & 0\\0 & 4\end{bmatrix}\longrightarrow\begin{bmatrix}4 & 0\\0 & 24\end{bmatrix}.
$$
The idea is again to use the Eulidean algorithm. After adding column 1 to column 2, you may have to do several row operations to obtain the GCD of the two original diagonal elements. In the example above, only one row operation was needed for this.
Addendum: Details for the particular example in the OP's question.
If you add row 2 to row 1, and then multiply row 1 by $-1$, you get a pivot of 1. You can then subtract suitable multiples of row 1 from rows 2 and 4, and end up with
$$
\begin{bmatrix}
1 & 561 & -174 & -80 \\
0 & -3477 & 1080 & 474 \\
0 & -255 & 81 & 24 \\
0 & -4182 & 1299 & 570
\end{bmatrix},
$$
which, by subtracting multiples of column 1 from columns 2, 3, and 4, leads to
$$
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & -3477 & 1080 & 474 \\
0 & -255 & 81 & 24 \\
0 & -4182 & 1299 & 570
\end{bmatrix}.
$$
We next need to do a series of row operations involving rows 2, 3, and 4 that results in the GCD of 3477, 255, and 4182. If we always choose the pivot of smallest nonzero magnitude, the steps would be
- Subtract 14 times row 3 from row 2, and 17 times row 3 from row 4.
- Add 2 times row 2 to row 3, and subtract row 2 from row 4.
- Subtract row 4 from row 2, and add row 4 to row 3.
- Add 3 times row 3 to row 2, and 6 times row 3 to row 4.
- Add row 2 to row 3, and subtract row 2 from row 4.
- Add 2 times row 3 to row 2.
- Swap rows 2 and 3.
- Negate row 2.
You should now have
$$
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 3 & 234 & -1410 \\
0 & 0 & -651 & 3906 \\
0 & 0 & -147 & 882
\end{bmatrix}.
$$
You can now subtract suitable multiples of column 2 from columns 3 and 4. Then you just have to deal with the $2\times2$ block in the lower right. You should be able to get it from here.
Best Answer
I'll answer your second question first.
Let
$$P(x) = \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&(x-2)^2&0\\0&0&0&(x-1)(x-2)\end{bmatrix}.$$
We need to fix the bottom right $2 \times 2$ principal submatrix. I have explained how to do that here, so I'll just present the results here:
\begin{align*} \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&(x-2)^2&0\\0&0&0&(x-1)(x-2)\end{bmatrix} &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&(x-2)^2&(x-1)(x-2)\\0&0&0&(x-1)(x-2)\end{bmatrix} \\ &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&(x-2)^2&(x-2)^2+(x-2)\\0&0&0&(x-1)(x-2)\end{bmatrix} \\ &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&(x-2)^2&x-2\\0&0&0&(x-1)(x-2)\end{bmatrix} \\ &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&x-2&(x-2)^2\\0&0&(x-1)(x-2)&0\end{bmatrix} \\ &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&x-2&0\\0&0&(x-1)(x-2)&-(x-1)(x-2)^2\end{bmatrix} \\ &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&x-2&0\\0&0&0&-(x-1)(x-2)^2\end{bmatrix} \\ &\mapsto \begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&x-2&0\\0&0&0&(x-1)(x-2)^2\end{bmatrix} \end{align*}
Irreducible factors are $x-2$ and $x-1$, while the elementary divisors are $x-2$, $(x-2)^2$ and $x-1$.
I think, but I'm not sure, that what you suggest in your first question can always be done to obtain the correct result. One would have to prove it to actually use it. Personally, I find it easier to just properly reduce the matrix polynomial or to fix it to become the Smith normal form, like I did above.