Roots and squarefree factorisation

ca_poly_roots computes the roots of a polynomial with exact coefficients, as exact numbers, together with their multiplicities. This is possible only when the roots can be expressed in the forms that ca_t supports, and the function says when it fails rather than returning approximations. ca_mat_eigenvalues is built on it.

Method#

The polynomial f must have a leading coefficient that is proved nonzero. ca_poly_roots first computes the squarefree factorisation f = c f_1^{e_1} \cdots f_k^{e_k} with ca_poly_factor_squarefree, which implements Yun's algorithm [Yun1976]: with g = \gcd(f, f'), the polynomials v = f/g and w = f'/g are updated by d = \gcd(v, w - v'), v \gets v/d, w \gets (w - v')/d, and the successive d are the squarefree factors of multiplicity 1, 2, \ldots. Every GCD needs zero tests of leading coefficients (see ca_poly_gcd), so this step can fail. Each squarefree factor is then solved:

• after removing the factor x^m if the constant terms are proved zero, a factor of degree 1 gives its root by a division;
• a factor whose coefficients are all rational is solved exactly over \overline{\mathbb{Q}} with qqbar_roots_fmpz_poly, which factors it over \mathbb{Z} and isolates the roots of each irreducible factor (see Roots of polynomials and eigenvalues); the roots become ca_t numbers with ca_set_qqbar, so they are elements of number fields generated by them;
• other factors of degree 2 and 3 are solved by the quadratic and cubic formulas, which introduce square roots, cube roots and the root of unity e^{2 \pi i/3};
• other factors of degree 4 or higher make the computation fail.

The roots are returned in no particular order. To build a polynomial from given roots, use ca_poly_set_roots, which multiplies the linear factors by binary splitting.

Example#

The roots of (x^2 + 1)^2 (x - \pi) are found with their multiplicities. The quintic x^5 - x - 1, which is not solvable by radicals, still has exact roots since its coefficients are rational: they are algebraic numbers of degree 5. The polynomial x^4 - \pi cannot be solved.

#include <flint/ca_poly.h>

int main(void)
{
    ca_ctx_t ctx;
    ca_poly_t f, g;
    ca_vec_t roots;
    ca_t c;
    ulong exp[5];
    slong i;

    ca_ctx_init(ctx);
    ca_poly_init(f, ctx);
    ca_poly_init(g, ctx);
    ca_vec_init(roots, 0, ctx);
    ca_init(c, ctx);

    /* f = (x^2 + 1)^2 (x - pi) */
    ca_poly_x(f, ctx);
    ca_poly_mul(f, f, f, ctx);
    ca_one(c, ctx);
    ca_poly_set_coeff_ca(f, 0, c, ctx);
    ca_poly_mul(f, f, f, ctx);
    ca_poly_x(g, ctx);
    ca_pi(c, ctx);
    ca_neg(c, c, ctx);
    ca_poly_set_coeff_ca(g, 0, c, ctx);
    ca_poly_mul(f, f, g, ctx);

    if (ca_poly_roots(roots, exp, f, ctx))
    {
        for (i = 0; i < ca_vec_length(roots, ctx); i++)
        {
            flint_printf("multiplicity %wu: ", exp[i]);
            ca_print(ca_vec_entry(roots, i), ctx);
            flint_printf("\n");
        }
    }

    /* x^5 - x - 1 has a nonsolvable Galois group, but rational coefficients */
    ca_poly_zero(g, ctx);
    ca_one(c, ctx);
    ca_poly_set_coeff_ca(g, 5, c, ctx);
    ca_neg_one(c, ctx);
    ca_poly_set_coeff_ca(g, 1, c, ctx);
    ca_poly_set_coeff_ca(g, 0, c, ctx);
    flint_printf("x^5 - x - 1: %d\n", ca_poly_roots(roots, exp, g, ctx));
    ca_printn(ca_vec_entry(roots, 0), 10, ctx); flint_printf("\n");

    /* x^4 - pi has no exact roots that ca_poly_roots can find */
    ca_pi(c, ctx);
    ca_neg(c, c, ctx);
    ca_poly_zero(g, ctx);
    ca_poly_set_coeff_ca(g, 0, c, ctx);
    ca_one(c, ctx);
    ca_poly_set_coeff_ca(g, 4, c, ctx);
    flint_printf("x^4 - pi: %d\n", ca_poly_roots(roots, exp, g, ctx));

    ca_poly_clear(f, ctx);
    ca_poly_clear(g, ctx);
    ca_vec_clear(roots, ctx);
    ca_clear(c, ctx);
    ca_ctx_clear(ctx);
    return 0;
}

multiplicity 1: 3.14159 {a where a = 3.14159 [Pi]}
multiplicity 2: 1.00000*I {a where a = I [a^2+1=0]}
multiplicity 2: -1.00000*I {-a where a = I [a^2+1=0]}
x^5 - x - 1: 1
1.167303978
x^4 - pi: 0