Multiplicar las raíces de polinomios utilizando sólo sus coeficientes enteros

Question

Estoy tratando de escribir una función que toma los coeficientes enteros de dos polinomios y devuelve los coeficientes de un polinomio que tiene una raíz para cada forma en que se puede multiplicar un par de raíces del par de polinomios de entrada.

La aplicación que tenía en mente es la representación exacta de números algebraicos. En lugar de almacenar un número dado $x$ como un flotador o un símbolo, almacenaría los coeficientes de un polinomio que tiene $x$ como raíz. Operar con números almacenados de esa manera, sin perder precisión, requiere operar con las raíces sin calcularlas realmente.

(Tengo claro que esta estrategia de representación tiene graves defectos, tanto en términos de eficiencia como de multiplicidad de soluciones, pero sigo queriendo resolver el problema).

Ejemplo

Supongamos que las listas de coeficientes de los polinomios de entrada son [1,0,-1] y [1,5,6]. Es decir, la entrada es el polinomio $x^2 - 1$ y el polinomio $x^2 + 5x + 6$ . Tienen raíces $-1,+1$ y $-2,-3$ respectivamente. Por lo tanto, el polinomio de salida debe tener raíces $-1 \cdot -2,-1 \cdot -3,1 \cdot -2,1 \cdot -3 = -3,-2,2,3$ y es igual a $(x-3)(x-2)(x+2)(x+3)$ que es $x^4 - 13x^2 + 36$ y así la salida es [1,0,-13,0,36].

Parcialmente resuelto

He trabajado en el problema, y he conseguido factorizar los efectos de los polinomios de entrada en piezas separadas (eso creo, al menos). Me llevaría demasiado tiempo explicar cómo llegué a donde estoy... basta con decir que fue muy ad-hoc.

La siguiente función parece calcular lo que quiero, y es casi independiente de las raíces:

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

Lo siento si no ha quedado del todo claro. He reducido el código todo lo posible.

He aquí una formulación alternativa de F, como función recursiva:

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

Como puedes ver, la única dependencia de las raíces es para calcular la función F. Parece que siempre se puede reducir a una expresión en términos de los coeficientes... Puedo hacerlo a mano para polinomios de grado pequeño, y puedo tropezar para reducirlo para los más grandes, pero no veo el patrón todavía...

Question

Así que...

• ¿Está ya resuelto este problema?
• ¿Hay alguna técnica que deba conocer y que pueda ayudar?
• ¿Me estoy perdiendo algo obvio?

Actualización

He utilizado la solución de Jyrki y la he puesto en práctica. Puedes encontrar el código en mi Repo de github de PolyNumber .

Answer 1

El problema se puede resolver en general utilizando la teoría clásica de resultantes . Para formar el polinomio cuyas raíces son productos $z=xy$ de raíces de $P(x)=0$ y $Q(y)=0$ , toma la resultante de $P(x)$ y $x^d Q(z/x)$ con respecto a $x$ , donde $d = \deg Q$ . En gp :

M(P,Q) = polresultant(P, x^poldegree(Q)*subst(Q,y,z/x), x)

y luego para su ejemplo

M(x^2-1, y^2+5*y+6)

devuelve $z^4 - 13z^2 + 36$ como debería.

Answer 2

Se sugiere el siguiente algoritmo. Sea $\alpha$ sea una raíz desconocida de un polinomio $f(x)$ de grado $n$ y $\beta$ sea una raíz desconocida de un polinomio $g(x)$ de grado $m$ .

Utilizando la ecuación $f(\alpha)=0$ nos permite reescribir cualquier poder $\alpha^k$ como una combinación lineal de $1,\alpha,\alpha^2,\ldots,\alpha^{n-1}$ . Del mismo modo, utilizando la ecuación $g(\beta)=0$ nos permite reescribir cualquier poder $\beta^\ell$ como un lineal combinación lineal de $1,\beta,\beta^2,\ldots,\beta^{m-1}$ . Uniendo estas dos piezas se demuestra que podemos escribir cualquier monomio $\alpha^k\beta^\ell$ como una combinación lineal del $mn$ cantidades $\alpha^i\beta^j, 0\le i<n, 0\le j<m.$ Denote $c_k=\alpha^i\beta^j$ , donde $k=mi+j$ oscila entre $0$ a $mn-1$ .

A continuación vamos a utilizar las expansiones de $(\alpha\beta)^t$ , $0\le t\le mn$ en términos de $c_k$ :s. Que estos sean $$ (\alpha\beta)^t=\sum_{k=0}^{mn-1}a_{kt}c_k. $$ Aquí los coeficientes $a_{tk}$ son números enteros. Buscamos encontrar enteros $x_t,t=0,1,\ldots,mn$ , de tal manera que $$ \sum_{t=0}^{mn}x_t(\alpha\beta)^t=0.\qquad(*) $$ Sustituyamos nuestra fórmula para la potencia $(\alpha\beta)^t$ arriba. La ecuación $(*)$ se convierte en $$ 0=\sum_t\sum_kx_ta_{tk}c_k=\sum_k\left(\sum_t x_ta_{kt}\right)c_k. $$ Esto será trivialmente cierto, si los coeficientes de todos los $c_k$ :s se desvanecen, es decir, es la ecuación $$ \sum_{t=0}^{mn}a_{kt}x_t=0 \qquad(**) $$ es válida para todos los $k=0,1,2,\ldots,mn-1$ . Aquí hay $mn$ ecuaciones lineales homogéneas en el $mn+1$ incógnitas $x_t$ . Por tanto, el álgebra lineal dice que tenemos garantizado el éxito en el sentido de que existe un vector no trivial $(x_0,x_1,\ldots,x_{mn})$ de números racionales que es una solución de $(**)$ . Además, multiplicando por el mínimo común múltiplo de los denominadores podemos hacer que todos los $x_t$ :s enteros.

El polinomio $$ F(x)=\sum_{t=0}^{mn}x_tx^t $$ es entonces una respuesta.

---

Hagamos su ejemplo de $f(x)=x^2-1$ y $g(x)=x^2+3x+2$ . Aquí $f(\alpha)=0$ nos dice que $\alpha^2=1$ . Del mismo modo, $g(\beta)=0$ nos dice que $\beta^2=-3\beta-2$ . Esto es todo lo que necesitamos de los polinomios. Aquí $c_0=1$ , $c_1=\beta$ , $c_2=\alpha$ y $c_3=\alpha\beta$ . Escribamos la potencia $(\alpha\beta)^t$ , $t=0,1,2,3,4$ en términos de $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. $$ Así pues, buscamos soluciones del sistema lineal homogéneo $$ \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. $$ Busquemos una solución con $x_4=1$ (si los polinomios con los que comenzó son mónicos, entonces esto siempre funciona AFAICT). La tercera ecuación nos dice que deberíamos tener $x_3=0$ . La segunda ecuación nos permite resolver que $x_2=-5$ . La última ecuación y nuestro conocimiento de $x_3$ dice que $x_1=0$ . La primera ecuación dice entonces que $x_0=4.$ Esto nos da la salida $$ x^4-5x^2+4. $$ Las posibilidades de $\alpha$ son $\pm1$ y las posibilidades de $\beta$ son $-1$ y $-2$ . Podemos comprobar fácilmente que los cuatro productos posibles $\pm1$ y $\pm2$ son ceros de este polinomio cuártico.

Mi solución al sistema lineal fue ad hoc, pero hay algoritmos conocidos para ello: eliminación, solución explícita en términos de menores,...

Answer 3

El problema se ha resuelto, entran en el polinomio simétrico:

http://en.wikipedia.org/wiki/Symmetric_polynomial#Relation_with_the_roots_of_a_monic_univariate_polynomial

Answer 4

Auto-respuesta. El código relevante de mi implementación actual, que sólo maneja polinomios enteros mónicos (sin incluir las toneladas de cosas de tipo matriz/polinomio):

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