Compute adjugate matrix over commutative ring Let $A$ be a $n\times n$ matrix over a commutative ring. I'm looking for a good method to compute its adjugate matrix.
My current approach is to use the Cayley-Hamilton theorem:
$$\text{adj}(A) = -(A^{n-1} + c_1A^{n-2}+\ldots +c_{n-2}A+c_{n-1}\text{I})
$$
where $$\lambda^n + c_1\lambda^{n-1} + \ldots + c_{n-1}\lambda +c_n = \det(\lambda\text{I}-A).$$
i.e. the $c_i$ are the coefficients of the characteristic polynomial, which I can determine using the Berkowitz algorithm. In particular, I never have to make divisions, so all works over an arbitrary commutative ring.
However, this method needs $\mathcal{O}(n^4)$ arithmetical ring operations, so there might be better alternatives. I hope for something like $\mathcal{O}(n^3)$, as for matrix inversion over a field.
 A: A reference dated 2001-2003 is http://www4.ncsu.edu/~kaltofen/bibliography/01/KaVi01.pdf
Over a commutative ring, the complexities of the calculation of $\det(A)$ and $adj(A)$ are the same. With the standard multiplication, the authors obtain a bit-complexity in $n^{10/3+\epsilon}N^{1+\epsilon}$ where $N$ is the length of the entries of $A$. Yet, the programming seems complicated ; there is a simpler method with exponent $3.5$.
Conclusion: If you have a robust and simple method with exponent $4$, then it is perfect.
A: (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.
