Hi "DC".  I think that I have worked out that the Hermite normal form is a "trichotomous quasipolynomial" in the variable $p$.  If $f:\mathbb{Z} \to \mathbb{Z}$ is a function, then my definition is that $f$ is a trichotomous quasipolynomial if it is a quasipolynomial for $x \gg 0$, possibly a different quasipolynomial for $x \ll 0$, and in between finitely many unrestricted values.

I think that if $R$ is a canonically Euclidean ring — Euclidean with canonically chosen quotients and remainders — then there is a Hermite normal form for matrices over $R$.  In particular, I think that $R$ does not have to be a Euclidean domain.

As a first try, let $A$ be the ring of all functions $f:\mathbb{Z} \to \mathbb{Z}$, using pointwise quotients and remainders.  Hermite normal form over this ring is a model of computing Hermite normal form for any $\mathbb{Z}$ family of integer matrices.   $A$ is sort-of a Euclidean ring, except it isn't Noetherian.

Let $B$ be the subring of $A$ consisting of trichotomous quasipolynomials.  <s>Then I believe that $B$ is Noetherian and it is a Euclidean ring.  If that is correct, then you obtain a Hermite normal form that is also a trichotomous quasipolynomial.</s>

----------

It's not correct, at least not in any obvious way.  It is easy to check that $B$ is not only a subring, but is also closed under quotients with pointwise remainders.  In order to compute how $b(x)$ divides into $a(x)$, you can reduce to the case in which $a(x)$ and $b(x)$ are both polynomials.  Then for starters there is a quotient and a remainder in $\mathbb{Q}[x]$:
$$a(x) = q(x)b(x) + r(x).$$
Since $r(x)$ has lower degree than $b(x)$, its values are smaller than those of $b(x)$ when $x \gg 0$, so that part is okay; but it may be negative and $q(x)$ may not be integral.  We can fix all that by rounding $q(x)$, which creates quasipolynomial behavior; and by changing it by 1 to make $r(x)$ positive.  Also since $r(x)$ might have odd degree, there may be a different quasipolynomial solution when $x \ll 0$.

The part that is either not true or far from obvious is why $B$ is Euclidean.  My argument for that fell apart.  However, I can still show that the Euclidean algorithm for finitely many elements $a_1,\ldots,a_n$ of $B$ terminates in a finite number of steps, and consequently that the Hermite normal form stabilizes in a finite number of steps with trichotomous quasipolynomial entries.  The proof is a two-stage induction.  The outer stage is the sum of the degrees of $a_1,\ldots,a_n$.  If $\deg a_j < \deg a_k$ for some $j$ and $k$, then dividing $a_j$ into $a_k$ reduces the degree of $a_k$.  On the other hand, suppose that the degrees are all equal.  Then we can pass to a congruence class for the input $x$ in $\mathbb{Z}$ and apply a linear change of variables so that the leading coefficients are all integers.  Then (in the inner induction) the Euclidean algorithm on these polynomials, for $x \gg 0$, amounts to the Euclidean algorithm on their coefficients.  It is important, in this inner inductive part, to only change the variable $x$ once; don't worry if the lower-order coefficients are not integers.  Eventually a $0$ is produced and the degrees once again decrease.

Alas, this is a very informal writeup, but this time I think that it works.