FPL already has support for computing Hermite normal form for integer
matrices. Here we add support to computing Smith normal form.
This is a preparation for Barvinok's algorithm. Given a polyhedron $P =
\{ x | Ax + b = 0, Cx + d \leq 0 \}$, we must find a particular solution
$x_0$ of $Ax + b = 0$ in order to project lower-dimensional polyhedra
into full-dimensional ones. This requires the Smith normal form of the
integer matrix $A$.
The implementation here follows the algorithm in
[wikipedia](https://en.wikipedia.org/wiki/Smith_normal_form#Algorithm).