(EDIT: there is a flaw when using $\det A$ below, as pointed out by @NicolasMalebranche.)

This answer expands on my comment in the hope to convince you that Gaussian elimination can be run formally without the need to ever simplify ratios. Your observation that **this only applies to an integral domain** is correct, though.

Consider the $3\times 3$ case for simplicity. Let us call $p_1,p_2,p_3$ the pivots encountered during Gaussian elimination, and suppose for simplicity that no row exchanges are needed. Let us suppose at first that all the pivots are nonzero.

At the first step, we multiply our matrix $A$ on the left by
$$
L_1 =
\begin{bmatrix}
1 & 0 & 0\\
-a_{21} & p_1 & 0\\
-a_{31} & 0 & p_1
\end{bmatrix},
$$
which zeroes out the entries in position (2,1) and (3,1). Then we premultiply by a matrix with form
$$
L_2 =
\begin{bmatrix}
1 & 0 & 0\\
0 & 1 & 0\\
0 & * & p_2
\end{bmatrix},
$$
to zero out the (3,2) entry. Let us call the resulting upper triangular matrix
$$
U=
\begin{bmatrix}
p_1 & * & *\\
0 & p_2 & *\\
0 & 0 & p_3
\end{bmatrix}.
$$
We have written $L_2L_1A=U$, so $\operatorname{Adj} A = (p_1p_2p_3)A^{-1} = (p_1p_2p_3) U^{-1}L_2L_1$.

Now I'll show that the result of this computation can be determined by a sort of back-substitution. Let $b$ be a generic column of $L_2L_1$. We write $Ux=b$ as
$$
\begin{bmatrix}
p_1 & * & *\\
0 & p_2 & *\\
0 & 0 & p_3
\end{bmatrix}
\begin{bmatrix}
x_1\\x_2\\x_3
\end{bmatrix}
=
\begin{bmatrix}
b_1\\b_2\\b_3
\end{bmatrix}
$$
We aim to solve this linear system not for $x$ but for $p_1p_2p_3x$.

From the third row, we can determine $p_3x_3$ without divisions. Multiplying the second equation by $p_3$, we can get $p_2p_3x_2$ without divisions. Multiplying the first equation by $p_2p_3$, we can get $p_1p_2p_3x_1$ without divisions.

So we just need to multiply by the remaining factors to obtain the entries of $p_1p_2p_3x$.

For a generic $n\times n$ matrix, everything should be perfectly analogous. One needs to precompute the products $p_1p_2\dots p_i$ and $p_i p_{i+1}\dots p_n$ for each $i$; once you do that, I think that all the computations can be performed in $O(n^3)$.

What if some pivots are zero? If we use total pivoting, the only breakdown case is when we arrive at a point when the last $n-k$ rows of $L_{k}\dots L_2L_1U$ are zero. We can restrict to the case $k=n-1$, otherwise $A$ has rank $n-2$ or lower and $\operatorname{Adj} A$ is trivially zero. So $p_n=0$ and all the previous pivots are nonzero. Set $L_n=I$, and factor as above to get $LA=U$, with $L$ lower triangular and invertible. We have $\operatorname{Adj} A = \operatorname{Adj} U \operatorname{Adj} L^{-1}$. The term $\operatorname{Adj} L^{-1}$ is a multiple of $L$, so it should pose little problems. The term $\operatorname{Adj} U$ has nonzero only in its last row (since the last row of $U$ is zeros), and it should be possible to compute these determinants in $O(n^2)$ each, since each is the determinant of an invertible triangular matrix + a rank-1 matrix (failing everything, use the matrix determinant lemma and the trick explained previously to compute $\det B B^{-1}u$, where $B$ is the top $(n-1)\times(n-1)$ block of $U$).

All of this is absolutely inelegant, I concur.

1more comment