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.

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.