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.



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....


