Multiplying the roots of polynomials using only their integer coefficients I am trying to write a function that takes the integer coefficients of two polynomials and returns the coefficients of a polynomial that has a root for each way you can multiply a root pair from the pair of input polynomial.
The application I had in mind is representing algebraic numbers exactly. Instead of storing a given number $x$ as a float or a symbol, I'd store the coefficients of a polynomial that has $x$ as a root. Operating on numbers stored that way, without losing precision, requires operating on the roots without actually computing the roots.
(It's clear to me that this representation strategy has serious flaws, both in terms of efficiency and multiplicity of solutions, but I still want to solve the problem.)
Example
Suppose the coefficient lists of the input polynomials are [1,0,-1] and [1,5,6]. That is to say, the input is the polynomial $x^2 - 1$ and the polynomial $x^2 + 5x + 6$. They have roots $-1,+1$ and $-2,-3$ respectively. Thus the output polynomial must have roots $-1 \cdot -2,-1 \cdot  -3,1 \cdot  -2,1 \cdot -3 = -3,-2,2,3$ and is equal to $(x-3)(x-2)(x+2)(x+3)$ which is $x^4 - 13x^2 + 36$ and so the output is [1,0,-13,0,36].
Partially Solved
I have worked on the problem, and managed to factor the effects of the input polynomials into separate pieces (I think so, anyways). It would take far too long to explain how I got where I am... suffice it to say it was very ad-hoc.
The following function seems to compute what I want, and is almost independent of the roots:
var decreasings = DecreasingSequencesOfSize(
    length: poly1.Degree, 
    total: indexOfOutputCoefficient, 
    max: poly2.Degree);

var outputCoefficient =
    (from coefRepeatSequence in decreasings
     let coefs2Factor = 
        coefRepeatSequence
        .Select(e => coefs2[e])
        .Product()
     let coefs1Factor = F(coefRepeatSequence, roots1, coefs1)
     select coefs2Factor*coefs1Factor
    ).Sum();

private static BigInteger F(IReadOnlyList<int> d, IReadOnlyList<BigInteger> roots, IReadOnlyList<BigInteger> coefs) {
    return (from p in d.DistinctPermutations()
            select p.Select((e, i) => BigInteger.Pow(roots[i], e)).Product()
           ).Sum();
}

I'm sorry if that's not totally clear. I cut down the code as much as possible.
Here's an alternate formulation of F, as a recursive function:
private static BigInteger F2(IReadOnlyList<int> d, IReadOnlyList<BigInteger> roots) {
    var zeroes = d.Count(e => e == 0);
    if (zeroes == d.Count) return 1;

    var min = d.Min();
    if (min > 0) {
        var mecs = d.Where(e => e > 0).Select(e => e - min).ToArray();
        return BigInteger.Pow(roots.Product(), min) * F2(mecs, roots);
    }

    var decs = d.Where(e => e > 0).ToArray();
    return (from x in roots.Choose(d.Count - zeroes)
            select F2(decs, x)
            ).Sum();
}

As you can see, the only dependence on the roots is for computing the function F. It does seem that it can always be reduced to an expression in terms of the coefficients... I can do it by hand for small degree polynomials, and can stumble towards reducing it for larger ones, but I don't see the pattern yet...
Question
So...


*

*Is this problem already solved?

*Is there a technique I should know about that might help?

*Am I missing something obvious?


Update
I used Jyrki's solution, and implemented it. You can find the code at my PolyNumber github repo.
 A: The following algorithm suggests itself. Let $\alpha$ be an unknown root of a polynomial $f(x)$ of degree $n$, and $\beta$ be an unknown root of a polynomial $g(x)$ of degree $m$.
Using the equation $f(\alpha)=0$ allows us to rewrite any power $\alpha^k$ as
a linear combination of $1,\alpha,\alpha^2,\ldots,\alpha^{n-1}$. Similarly using the equation $g(\beta)=0$ allows us to rewrite any power $\beta^\ell$ as a linear
combination of $1,\beta,\beta^2,\ldots,\beta^{m-1}$. Putting these two pieces together shows that we can write any monomial $\alpha^k\beta^\ell$ as a linear combination of the $mn$ quantities $\alpha^i\beta^j, 0\le i<n, 0\le j<m.$
Denote $c_k=\alpha^i\beta^j$, where $k=mi+j$ ranges from $0$ to $mn-1$.
Next let us use the expansions of $(\alpha\beta)^t$, $0\le t\le mn$ in terms of the $c_k$:s. Let these be
$$
(\alpha\beta)^t=\sum_{k=0}^{mn-1}a_{kt}c_k.
$$
Here the coefficients $a_{tk}$ are integers. We seek to find integers $x_t,t=0,1,\ldots,mn$, such that
$$
\sum_{t=0}^{mn}x_t(\alpha\beta)^t=0.\qquad(*)
$$
Let us substitute our formula for the power $(\alpha\beta)^t$ above. The equation $(*)$ becomes
$$
0=\sum_t\sum_kx_ta_{tk}c_k=\sum_k\left(\sum_t x_ta_{kt}\right)c_k.
$$
This will be trivially true, if the coefficients of all $c_k$:s vanish, that is, is the equation
$$
\sum_{t=0}^{mn}a_{kt}x_t=0 \qquad(**)
$$
holds for all $k=0,1,2,\ldots,mn-1$. Here there are $mn$ linear homogeneous equations on the $mn+1$ unknowns $x_t$. Therefore linear algebra says that we are guaranteed to succeed in the sense that there exists a non-trivial vector
$(x_0,x_1,\ldots,x_{mn})$ of rational numbers that is a solution of $(**)$. Furthermore, by multiplying with the least common multiple of the denominators,
we can make  all the $x_t$:s integers.
The polynomial
$$
F(x)=\sum_{t=0}^{mn}x_tx^t
$$
is then an answer.

Let's do your example of $f(x)=x^2-1$ and $g(x)=x^2+3x+2$. Here $f(\alpha)=0$
tells us that $\alpha^2=1$. Similarly $g(\beta)=0$ tells us that $\beta^2=-3\beta-2$. This is all we need from the polynomials. Here $c_0=1$, $c_1=\beta$,
$c_2=\alpha$ and $c_3=\alpha\beta$. Let us write the power $(\alpha\beta)^t$,
$t=0,1,2,3,4$ in terms of the $c_k$:s.
$$
(\alpha\beta)^0=1=c_0.
$$
$$
(\alpha\beta)^1=\alpha\beta=c_3.
$$
$$
(\alpha\beta)^2=\alpha^2\beta^2=1\cdot(-3\beta-2)=-2c_0-3c_1.
$$
$$
(\alpha\beta)^3=\alpha\beta(-3\beta-2)=\alpha(-3\beta^2-2\beta)=\alpha(9\beta+6-2\beta)=\alpha(7\beta+6)=6c_2+7c_3.
$$
$$
(\alpha\beta)^4=(-3\beta-2)^2=9\beta^2+12\beta+4=(-27\beta-18)+12\beta+4=-14-15\beta=-14c_0-15c_1.
$$
We are thus looking for solutions of the homogeneous linear system
$$
\left(\begin{array}{ccccc}
1&0&-2&0&-14\\
0&0&-3&0&-15\\
0&0&0&6&0\\
0&1&0&7&0
\end{array}\right)
\left(\begin{array}{c}x_0\\x_1\\x_2\\x_3\\x_4\end{array}\right)=0.
$$
Let us look for a solution with $x_4=1$ (if the polynomials you started with
are monic, then this always works AFAICT). The third equation tells us that we should have $x_3=0$. The second equation allows us to solve that $x_2=-5$. The last equation and our knowledge of $x_3$ tells that $x_1=0$. The first equation
then tells that $x_0=4.$ This gives us the output
$$
x^4-5x^2+4.
$$
The possibilities for $\alpha$ are $\pm1$ and the possibilities for $\beta$ are $-1$ and $-2$. We can easily check that all the possible four products $\pm1$ and $\pm2$ are zeros of this quartic polynomial.
My solution to the linear system was ad hoc, but there are known algorithms for that: elimination, an explicit solution in terms of minors,...
A: The problem can be solved in general using the classical theory of
resultants.
To form the polynomial whose roots are products $z=xy$ 
of roots of $P(x)=0$ and $Q(y)=0$, take the resultant of
$P(x)$ and $x^d Q(z/x)$ with respect to $x$, where $d = \deg Q$.
In gp:
M(P,Q) = polresultant(P, x^poldegree(Q)*subst(Q,y,z/x), x)

and then for your example
M(x^2-1, y^2+5*y+6)

returns $z^4 - 13z^2 + 36$ as it should.
A: The problem has been solved, they fall under symmetric polynomial:
http://en.wikipedia.org/wiki/Symmetric_polynomial#Relation_with_the_roots_of_a_monic_univariate_polynomial
A: Self-answer. The relevant code of my current implementation, which only handles monic integer polynomials (not including the tons of matrix/polynomial type stuff):
private static IEnumerable<IntPolynomial<XTerm>> RewritePolyBasisUsing(XTerm termToRewrite,
                                                                       IntPolynomial<XTerm> rewrittenTerm) {
    var X = new XTerm(1);
    var r = rewrittenTerm - termToRewrite;

    var cur = (IntPolynomial<XTerm>)1;
    yield return cur;

    for (var i = 1; i < termToRewrite.XPower; i++) {
        cur *= X;
        yield return cur;
    }
    while (true) {
        cur = cur*X;
        var factor = cur.Coefficient(termToRewrite);
        cur += r*factor;
        yield return cur;
    }
}
private static IEnumerable<IntPolynomial<XTerm>> RewritePolyBasisUsing(IntPolynomial<XTerm> poly) {
    var termToRewrite = poly.Coefficients.MaxBy(e => e.Key.XPower);
    if (termToRewrite.Value != 1) throw new NotImplementedException();
    var rewrittenTerm = termToRewrite - poly;
    return RewritePolyBasisUsing(termToRewrite.Key, rewrittenTerm);
}

public static IntPolynomial<XTerm> MultiplyRoots(this IntPolynomial<XTerm> poly1,
                                                 IntPolynomial<XTerm> poly2) {
    var deg1 = poly1.Degree();
    var deg2 = poly2.Degree();

    var maxDegreeOfResult = deg1*deg2;
    if (maxDegreeOfResult == 0) return 1;

    var rewrittenPowers1 = RewritePolyBasisUsing(poly1).Select(p => p.ToPolynomialOverVariable1Of2());
    var rewrittenPowers2 = RewritePolyBasisUsing(poly2).Select(p => p.ToPolynomialOverVariable2Of2());
    var rewrittenPowerPairs =
        rewrittenPowers1.Zip(rewrittenPowers2,
                             IntPolynomial<XYTerm>.Multiply)
                        .Take(maxDegreeOfResult + 1);

    var linearSystem = IntMatrix.FromColumns(
        from col in rewrittenPowerPairs
        select from row in maxDegreeOfResult.Range()
               let exponentForX = row/deg2
               let exponentForY = row%deg2
               select col.Coefficient(new XYTerm(exponentForX, exponentForY)));

    var solvedSystem = linearSystem.Reduced();

    var degreeOfSolution = solvedSystem.Rows.Count(row => row.Any(cell => cell != 0));

    var lowCoefficients = solvedSystem.Columns[degreeOfSolution].Take(degreeOfSolution);

    return Polynomial.XToThe(degreeOfSolution) - lowCoefficients.TimesPolynomialBasis();
}

