How to obtain the rational solution of a linear system efficiently? Suppose that I have a linear system $AX=b$ with $A\in\mathbb{Z}^{n\times m}$ and $b\in\mathbb{Z}^{n}$. Assume that $AX=b$ has exactly one rational solution. Then how can I obtain this solution efficiently?
I know LU decomposition is a choice for this purpose. But it seems that the denominators may expand in the process of computing a LU decomposition, which I hope to avoid if the rational solution has small denominators.
If I already have a numerical solution of $AX=b$, does it help?
 A: The following answer is only theoretical. Suppose you first reduce your matrix to Smith normal form $D = S A T$, where $S$ and $T$ are integral matrices with integral inverses (the determinants of $S$ and $T$ are each either $+1$ or $-1$) and $D$ is diagonal, whose diagonal elements are the so-called elementary divisors of $A$. Then you reduce your problem to solving $DY=c$, where $c = S b$ and $X = T Y$. Since by $c$ is still integral and $D$ is now diagonal, the denominators in $Y$ will be no larger than the largest entry in $D$ (the largest elementary divisor of $A$). Then, since $X$ is obtained from $Y$ by multiplication with another integral matrix, the denominators in $X$ will be no larger than the lcm of the denominators in $Y$.
Unfortunately, I don't think that reduction to Smith normal form is a cheap operation, so whether it is worth going through that step in your problem will depend on other factors.
EDIT: Actually, the Smith normal form is overkill for your problem. The Hermite normal form is basically the integer version of the LU decomposition that you already had in mind: $U = S A$, where $U$ is integral upper triangular, and $S$ is integral with integral inverse. The new system you need to solve is $U X = S b$, which can be done by back substitution, where at each step you only divide by the corresponding diagonal element of $U$. Hence, the final denominators in $X$ will be no larger than the lcm of the diagonal elements of $U$.
The linked Wikipedia article on the Hermite normal form mentions that modern algorithms for its computation have polynomially bounded time and space complexity (respectively, in terms of the matrix size and logarithmic size of the matrix entries). If you (or your computer algebra system) use such an algorithm, it might be your best bet for provably avoiding intermediate expression swell. The practicality of this method for your problem of course would still be subject to experimentation.
EDIT 2: I've also found a very well written survey of several algorithms for solving rational linear systems, with practical timing comparisons in the case of sparse systems.

Cook, William; Steffy, Daniel E., Solving very sparse rational systems of equations, ACM Trans. Math. Softw. 37, No. 4, Paper No. 39, 21 p. (2011). ZBL1365.65121. (author's copy)

A: This is standard and can be done using $p$-adic lifting, see for example the following article:

*

*J.D. Dixon, Exact solution of linear equations using $P$-adic expansions, Numer. Math. 40 (1982) pp. 137–141, doi:10.1007/BF01459082.

This also exploits that there is a unique solution.
