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. 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.
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.
