Computing determinants of matrices of linear forms
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....
(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.
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.