Roots modulo a prime and modulo composite numbers

nmod_poly_roots finds all roots of a nonzero polynomial f in the prime field \mathbb{F}_p and returns them as the linear factors x - a of an nmod_poly_factor_t, with the multiplicity of each root as its exponent if requested. nmod_poly_roots_factored does the same over \mathbb{Z}/n\mathbb{Z} for a composite modulus n whose prime factorisation is known. Over \mathbb{Z}/n\mathbb{Z} the number of roots is not bounded by the degree, and the multiplicity of a is defined as the largest e with (x - a)^e \mid f in (\mathbb{Z}/n\mathbb{Z})[x].

To check whether a polynomial splits into distinct linear factors with nonzero roots, and find them if it does, the specialised nmod_poly_find_distinct_nonzero_roots in nmod_poly is faster, since it skips the squarefree factorisation (see Finding roots that are expected to be distinct).

Algorithms#

Modulo a prime, nmod_poly_roots computes the squarefree factorisation of f, then for each squarefree part g the two gcds \gcd(g, x^{(p-1)/2} - 1) and \gcd(g, x^{(p-1)/2} + 1), whose product is the product of the distinct nonzero linear factors of g, and splits them with random shifts \gcd(u, (x + c)^{(p-1)/2} - 1). The cost is dominated by the exponentiations, O(M(n) \log p \log n) in expectation for degree n; the irreducible factors of degree at least 2 are never separated from each other. The root 0 is split off directly.

Modulo n = \prod_j p_j^{k_j}, nmod_poly_roots_factored finds the roots modulo each p_j, lifts every one of them to all roots modulo p_j^{k_j} that reduce to it, using Lemma 3 of Roots of polynomials over finite fields with doubling precision, and combines the roots for different primes with the Chinese remainder theorem. A root that is simple modulo p_j lifts uniquely, while a multiple root can lift to many roots or to none. The function returns 0, without roots, if a list of roots would exceed 2^{32} entries on 64-bit systems.

Example#

Modulo 7, f = (x - 1)^3 (x - 2)(x^2 + 1) has the roots 1 and 2 (the factor x^2 + 1 is irreducible since 7 \equiv 3 \bmod 4). Modulo 120 = 2^3 \cdot 3 \cdot 5, x^2 - 1 has 4 \cdot 2 \cdot 2 = 16 roots, and x^2 has the roots a with a^2 \equiv 0, that is, the multiples of 2^2 \cdot 3 \cdot 5 = 60:

#include <stdlib.h>
#include <flint/flint.h>
#include <flint/ulong_extras.h>
#include <flint/nmod.h>
#include <flint/nmod_poly.h>
#include <flint/nmod_poly_factor.h>

static int cmp(const void * a, const void * b)
{
    ulong x = *(const ulong *) a, y = *(const ulong *) b;
    return (x > y) - (x < y);
}

/* print the roots a (from the factors x - a) in increasing order,
   with the multiplicity if it is greater than 1 */
static void print_roots(const nmod_poly_factor_t r, ulong n)
{
    ulong a[64][2];
    slong i;
    for (i = 0; i < r->num; i++)
    {
        a[i][0] = (n - r->p[i].coeffs[0]) % n;
        a[i][1] = r->exp[i];
    }
    qsort(a, r->num, sizeof(a[0]), cmp);
    flint_printf("%wd roots mod %wu:", r->num, n);
    for (i = 0; i < r->num; i++)
    {
        flint_printf(" %wu", a[i][0]);
        if (a[i][1] > 1)
            flint_printf("^%wu", a[i][1]);
    }
    flint_printf("\n");
}

int main(void)
{
    nmod_poly_t f, g;
    nmod_poly_factor_t r;
    n_factor_t fac;
    slong i;

    nmod_poly_init(f, 7);
    nmod_poly_init(g, 7);
    nmod_poly_factor_init(r);

    nmod_poly_set_str(f, "2 7  6 1");          /* x - 1 */
    nmod_poly_pow(f, f, 3);
    nmod_poly_set_str(g, "2 7  5 1");          /* x - 2 */
    nmod_poly_mul(f, f, g);
    nmod_poly_set_str(g, "3 7  1 0 1");        /* x^2 + 1 */
    nmod_poly_mul(f, f, g);

    nmod_poly_roots(r, f, 1);
    for (i = 0; i < r->num; i++)
        flint_printf("root %wu with multiplicity %wd\n",
            nmod_neg(r->p[i].coeffs[0], f->mod), r->exp[i]);

    nmod_poly_clear(f);
    nmod_poly_init(f, 120);
    n_factor_init(&fac);
    n_factor(&fac, 120, 1);

    nmod_poly_set_str(f, "3 120  119 0 1");    /* x^2 - 1 */
    if (nmod_poly_roots_factored(r, f, 1, &fac))
        print_roots(r, 120);

    nmod_poly_set_str(f, "3 120  0 0 1");      /* x^2 */
    if (nmod_poly_roots_factored(r, f, 1, &fac))
        print_roots(r, 120);

    nmod_poly_clear(f);
    nmod_poly_clear(g);
    nmod_poly_factor_clear(r);
    return 0;
}

root 2 with multiplicity 1
root 1 with multiplicity 3
16 roots mod 120: 1 11 19 29 31 41 49 59 61 71 79 89 91 101 109 119
2 roots mod 120: 0^2 60^2

The example reads only the constant coefficients of the factors, which have the modulus of f. Both roots of x^2 have multiplicity 2, since (x - 60)^2 = x^2 - 120x + 3600 \equiv x^2 \pmod{120}.