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.

up vote 5 down vote accepted

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

  • Well, that was fast. I'll try this out and report back. – Craig Gidney Sep 10 '13 at 17:11
  • 1
    Yup, finished implementing it and it passes tests. Will link to the repo when that's up.. – Craig Gidney Sep 11 '13 at 17:41
  • @Strilanc: Glad to hear that you succeeded in implementing this yourself. – Jyrki Lahtonen Sep 11 '13 at 17:51
  • I should note that I did run into an issue with multiplicity of roots. In some cases a double-root would reduce the system by one degree of freedom, and in other cases it wouldn't. I ended up having to test against upper and lower bounds on 'divides' instead of exact match. – Craig Gidney Sep 11 '13 at 18:14
  • 1
    The tricky part is that it can happen when they have no roots in common. For example: (x-2)(x-6) vs (x-3)(x-15) creates (x-6)(x-18)(x-30)^2. I think it being created by the combination has to do with when it is and isn't removing a degree of freedom. For my application every redundant root lost is actually a benefit so I'm not worried about it. – Craig Gidney Sep 11 '13 at 20:09

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.

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

  • Could you expand on the solution? The section of the article you linked only seems to cover facts I was already using anyways, like how the i'th coefficient is the sum of all ways you can choose i roots to multiply together or that a solution should exist. I'll read some more. – Craig Gidney Sep 10 '13 at 17:01
  • If you have the coefficients, like you do, you can use the discrete convolution operator, it is dead simple to implement. Given that you have roots 1,2,3,4,5,6 you have the coefficients [1 -1],[1 -2],[1 -3],[1 -4],[1 -5],[1 -6], respectively. Now convolve (((([1 -1]*[1 -2])*[1 -3])*[1 -4])*[1-5])*[1 -6]. This does not require the polynomials to be monic or even order 1. – riboch Sep 10 '13 at 21:47

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();
}

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.