Computing determinants of matrices of linear forms Suppose we have three $n \times n$ matrices $A$, $B$, $C$ with floating point entries. We would like to compute the polynomial $\det (xA+yB+zC)$. At least in Mathematica, and I think in all computer algebra systems, this will take $n!$ steps; Mathematica chokes around $n=15$. 
What's the smart way to do this? I have some ideas, but they are all complicated enough that I don't want to implement them before hearing from others.
This question is on the border between MO, cstheory and Mathematica, but I suspect that this is the right place to start.
ADDED IN RESPONSE TO QUESTIONS BELOW The trouble with Gaussian elimination is that you have to divide by polynomials in $(x,y,z)$, and the expressions soon get huge. Interpolation, probably at roots of unity so that the interpolation matrix will be unitary, is my best idea, but I wanted to see if there was a better one before I tried it.
 A: The smart way of doing this is by interpolation (this, by the way, is how Mathematica computes the characteristic polynomial), and in Mathematica you can use InterpolatingPolynomial[]  to do this in a couple-of-line program.
Remark In fact, to compute the characteristic polynomial of an integer matrix, Mathematica uses interpolation twice, in effect. That's because to compute the ordinary determinant, it uses Chinese remaindering....
A: (1) Look up Gauss-Bareiss and Dodgson condensation - the determinant can be computed fraction-free over a PID.
(2) Try interpolation. I've used this in anger for the two-variable case, it does look
less pleasant with three variables.
A: To just add to the answers already here, it is possible to extract the individual coefficients of this polynomial exactly using the following method. Let
\begin{equation}
h(x) = \sum_{j = 0}^d c_j x^j
\end{equation}
be a polynomial of degree $d$, then
\begin{equation}
c_j = \frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} h(e^{\frac{2 \pi i k}{d + 1}})
\end{equation}
since
\begin{align}
\frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} h(e^{\frac{2 \pi i k}{d + 1}}) &= \frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} \left(\sum_{l = 0}^d c_l e^{\frac{2 \pi i k l}{d + 1}} \right) \\
&= \sum_{l = 0}^d \left(\frac{1}{d + 1} \sum_{k = 0}^d e^{\frac{2 \pi i (l - j) k}{d + 1}} \right) c_l \\
&= \sum_{l = 0}^d \delta_{jl} c_l \\
\frac{1}{d + 1} \sum_{k = 0}^d e^{-\frac{2 \pi i j k}{d + 1}} h(e^{\frac{2 \pi i k}{d + 1}}) &= c_j
\end{align}
where $\delta_{jl}$ is the Kronecker delta and $i$ is the imaginary unit. Letting 
\begin{equation}
f(x, y, z) = \det(x A + y B + z C)
\end{equation}
we can expand $f$ in powers of $x, y, z$ using multilinearity of the rows or columns of the matrix inside the determinant as
\begin{equation}
f(x, y, z) = \sum_{j = 0}^n x^j g_j(y, z) = \sum_{j = 0}^n x^j \left(\sum_{k = 0}^{n - j} a_{j, k} y^k z^{n - j - k}\right)
\end{equation}
Note that, since one choice of $x, y, z$ must be made for every row or column in the expansion, the powers of these variables must sum to $n$ in every term. This polynomial is completely specified by the coefficients $\{a_{j, k}\}$, and to extract them, it is sufficient to only evaluate $f$ at $z = 1$
\begin{equation}
f(x, y, 1) = \sum_{j = 0}^n x^j g_j(y, 1) = \sum_{j = 0}^n \sum_{k = 0}^{n - j} a_{j, k} x^j y^k
\end{equation}
One query of the polynomial 
\begin{equation}
g_j(y, 1) = \sum_{k = 0}^{n - j} a_{j, k} y^k
\end{equation}
is given by the above formula
\begin{equation}
g_j(y, 1) = \frac{1}{n + 1} \sum_{p = 0}^n e^{-\frac{2 \pi i j p}{n + 1}} f(e^{\frac{2 \pi i p}{n + 1}}, y, 1)
\end{equation}
Similarly,
\begin{equation}
a_{j,k} = \frac{1}{n - j + 1} \sum_{q = 0}^{n - j} e^{-\frac{2 \pi i k q}{n - j + 1}} g_j(e^{\frac{2 \pi i q}{n - j + 1}}, 1)
\end{equation}
Therefore
\begin{equation}
\det(x A + y B + z C) = \sum_{j = 0}^n \sum_{k = 0}^{n - j} a_{j,k} x^j y^k z^{n - j - k}
\end{equation}
where
\begin{equation}\
a_{j, k} = \frac{1}{(n + 1)(n - j + 1)}\sum_{p = 0}^{n} \sum_{q = 0}^{n - j} e^{-2 \pi i \left(\frac{jp}{n + 1} + \frac{kq}{n - j + 1} \right)} \det(e^{\frac{2 \pi i p}{n + 1}} A + e^{\frac{2 \pi i q}{n - j + 1}} B + C)
\end{equation}
requiring the evaluation of only $O(n^4)$ determinants to completely characterize the polynomial.
