It is "well known" that a matrix system $Ax=b$ where $A\in \Bbb Z^{m\times n}$, $x\in \Bbb Z^n,b\in\Bbb Z^m$ for some $m,n \in \Bbb N$, can be solved in polynomial time, using Smith/Hermite Normal Form computation. By "solved" I mean either "decide, for a given $A,b$, whether a solution exists" or "decide whether a solution exists and output one". My question is:

Can we claim anything more specific than

polynomial, such as cubic/quartic in $(m+n+\max{|A_{ij}}|)$? Or at least cubic/quartic if we assume that any arithmetic operation can be performed in unit time (ignoring possible blowup of the coefficients)?

The closest result I could find is in "Near optimal algorithms for computing Smith normal forms of integer matrices" by Storjohann where a nice complexity bound is presented for computing the Smith normal form, but it is not clear from the paper how to recover the unimodular transformation matrices (authors just write "In the future, we will present ... algorithm that compute (them)").

Another close result I found is this paper where a cubic algorithm is presented under the additional assumtion that we can use a unit-time oracle for computing greatest common divisor.

Another reference is this paper which deals with the case of non-singular square matrices. (Is there a simple reduction from a general $m\times n$ matrix to the non-singular square case?)

Another reference is this paper (thanks @Dima Pasechnik) which gives upper bounds on Hermite normal form computation, at least for the non-singular case: in the last section there is an outline how to generalize it for any matrices, but I again don't see how to recover the transformation matrix, as their algorithm does not use only elementary row/column operations.