Complexity of solving systems of linear diophantine equations 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.
 A: An improvement of Storjohann's result is due to 
Micciancio and Warinschi,
A linear space algorithm for computing the Hermite normal form. 
As to recovery of transformation matrices, in the classical algorithm one merely has to trace the elementary steps of the SNF/HNF computation; e.g. swapping rows/columns means that you have to multiply one of these matrices by the corresponding transposition; adding a row to another row means that you have to multiply one of them by the corresponding tranvection, etc. Note that these transformations are very sparse, i.e. these updates take much less time than general matrix multiplication.
Here it's less clear, and the computation is done modulo a number of primes, and this might give extra problems.
Storjohann and Labahn in Asymptotically fast computation of Hermite normal forms of integer matrices show how to compute unimodular $U$ so that $UA=H$, $H$ the HNF of $A$. There they require full column rank of $A$.
A: Ok, apparently, professor Storjohann sent me the reference I was looking for: Theorem 19 in A Fast Practical Deterministic Algorithm for Triangularizing Integer Matrices. (Thanks)
