Modular composition

Modular composition computes f(g) \bmod h for polynomials f, g, h with \deg f, \deg g < \deg h = n. It is the key operation of polynomial factorisation over finite fields, where it computes the iterates x^{p^i} \bmod h of the Frobenius map, and of the computation of minimal polynomials and of field embeddings. Horner's rule needs n - 1 multiplications modulo h, that is, O(n M(n)) operations; the algorithm of Brent and Kung [BrentKung1978] reduces this to O(\sqrt n M(n) + n^{(\omega + 1)/2}), where \omega is the exponent of matrix multiplication.

The Brent–Kung algorithm#

Let m = \lfloor \sqrt n \rfloor + 1 and write f = \sum_{j} f_j(x)\, x^{jm} with blocks f_j of length m. Then

The baby steps g^0, \ldots, g^{m-1} \bmod h form the rows of an m \times n matrix A, and all f_j(g) together are the rows of the product C = BA, where row j of the m \times m matrix B holds the coefficients of f_j. The giant steps evaluate \sum_j C_j H^j for H = g^m \bmod h by Horner's rule or rectangular splitting (_nmod_poly_mod_matrix_rows_evaluate), with about \sqrt n multiplications modulo h. The matrix product is done by nmod_mat_mul (see Multiplying matrices over Z/nZ).

The matrix A depends only on g and h. nmod_poly_precompute_matrix computes it once (it must be initialised with m rows and n columns), and nmod_poly_compose_mod_brent_kung_precomp_preinv then composes any f with g at the cost of the matrix product and the giant steps. The preinv versions take the inverse of the reversal of h, as explained in Computing modulo a fixed polynomial.

nmod_poly_compose_mod uses Horner's rule for \deg h < 7 and when \deg f \ge \deg h (Brent–Kung requires \deg f < \deg h), and the Brent–Kung algorithm otherwise.

Example: Rabin's irreducibility test#

A polynomial h of degree n over \mathbb{F}_p is irreducible if and only if x^{p^n} \equiv x \pmod h and \gcd(x^{p^{n/q}} - x, h) = 1 for every prime q dividing n [Rab1980b]. The powers x^{p^i} are computed by one exponentiation, x^p \bmod h, followed by repeated composition with x^p, using one precomputed matrix:

#include <flint/flint.h>
#include <flint/ulong_extras.h>
#include <flint/nmod_mat.h>
#include <flint/nmod_poly.h>
#include <flint/nmod_poly_factor.h>

/* x^(p^i) mod h for i = 0, ..., n, by composition with x^p */
int is_irreducible_rabin(const nmod_poly_t h)
{
    slong n = nmod_poly_degree(h), i, m = n_sqrt(n) + 1;
    ulong p = h->mod.n;
    nmod_poly_t hinv, xp, t, g;
    nmod_mat_t A;
    int result = 1;

    nmod_poly_init_mod(hinv, h->mod);
    nmod_poly_init_mod(xp, h->mod);
    nmod_poly_init_mod(t, h->mod);
    nmod_poly_init_mod(g, h->mod);
    nmod_mat_init(A, m, n, p);

    nmod_poly_reverse(hinv, h, h->length);
    nmod_poly_inv_series(hinv, hinv, h->length);

    nmod_poly_powmod_x_ui_preinv(xp, p, h, hinv);       /* x^p mod h */
    nmod_poly_precompute_matrix(A, xp, h, hinv);

    nmod_poly_set(t, xp);                               /* t = x^(p^1) */
    for (i = 1; i <= n; i++)
    {
        /* gcd test for the maximal proper divisors i = n/q of n */
        if (i < n && n % i == 0 && n_is_prime(n / i))
        {
            nmod_poly_set_coeff_ui(g, 1, 1);
            nmod_poly_sub(g, t, g);                     /* x^(p^i) - x */
            nmod_poly_gcd(g, g, h);
            if (!nmod_poly_is_one(g))
                result = 0;
            nmod_poly_zero(g);
        }
        if (i < n)                                      /* t = t(x^p) */
            nmod_poly_compose_mod_brent_kung_precomp_preinv(t, t, A, h, hinv);
    }
    if (!nmod_poly_is_gen(t))                           /* x^(p^n) == x */
        result = 0;

    nmod_mat_clear(A);
    nmod_poly_clear(hinv);
    nmod_poly_clear(xp);
    nmod_poly_clear(t);
    nmod_poly_clear(g);
    return result;
}

int main(void)
{
    flint_rand_t state;
    nmod_poly_t a, b, h;
    ulong p = 1000003;

    flint_rand_init(state);
    nmod_poly_init(a, p);
    nmod_poly_init(b, p);
    nmod_poly_init(h, p);

    nmod_poly_randtest_monic_irreducible(a, state, 61);   /* degree 60 */
    nmod_poly_randtest_monic_irreducible(b, state, 41);   /* degree 40 */

    flint_printf("a irreducible: %d (library: %d)\n",
                 is_irreducible_rabin(a), nmod_poly_is_irreducible(a));
    nmod_poly_mul(h, a, b);
    flint_printf("a*b irreducible: %d\n", is_irreducible_rabin(h));
    nmod_poly_mul(h, a, a);
    flint_printf("a^2 irreducible: %d\n", is_irreducible_rabin(h));

    nmod_poly_clear(a);
    nmod_poly_clear(b);
    nmod_poly_clear(h);
    flint_rand_clear(state);
    return 0;
}

a irreducible: 1 (library: 1)
a*b irreducible: 0
a^2 irreducible: 0

nmod_poly_is_irreducible is declared in nmod_poly_factor.h. Each composition costs one m \times m by m \times n matrix product and about m modular multiplications, against \log_2 p \approx 20 squarings modulo h for an exponentiation, which is where modular composition pays off for large p [GatSho1992] [KalSho1998].