Computing $J$-ideals of a matrix over a principal ideal domain

Given a square matrix $B$ over a principal ideal domain $D$ and an ideal $J$ of $D$, the $J$-ideal of $B$ consists of the polynomials $f\in D[X]$ such that all entries of $f(B)$ are in $J$. It has been shown that in order to determine all $J$-ideals of $B$ it suffices to compute a generating set of the $(p^t)$-ideal of $B$ for finitely many prime powers $p^t$. Moreover, it is known that $(p^t)$-ideals are generated by polynomials of the form $p^{t-s}\nu_s$ where $\nu_s$ is a monic polynomial of minimal degree in the $(p^s)$-ideal of $B$ for some $s\le t$. However, except for the case of diagonal matrices, it was not known how to determine these polynomials explicitly. We present an algorithm which allows to compute the polynomials $\nu_s$ for general square matrices. Exploiting one of McCoy's theorems we first compute some set of generators of the $(p^s)$-ideal of $B$ which then can be used to determine $\nu_s$. This algorithmic computation significantly extends our understanding of the $J$-ideals of $B$.