Computing modulo a fixed polynomial

Many algorithms compute in the quotient ring $(\mathbb{Z}/n\mathbb{Z})[x]/(f)$: arithmetic in a finite field $\mathbb{F}_{p^d} = \mathbb{F}_p[x]/(f)$ for an irreducible $f$ of degree $d$, distinct-degree factorisation, primality tests for polynomials, and the computation of $x^{p^i} \bmod f$ in general. An element of the quotient ring is represented by its remainder modulo $f$, a polynomial of degree less than $d = \deg f$, and the basic operation is a product followed by a reduction modulo $f$. The leading coefficient of $f$ must be invertible; for a field, $n$ must be prime and $f$ irreducible.

Reduction with a precomputed inverse#

A remainder modulo $f$ can be computed with nmod_poly_rem, which for large $d$ uses Newton iteration to find the quotient. That iteration computes an inverse of the reversed divisor each time. When the same $f$ is used many times, this inverse should be computed once. Write $\operatorname{rev}_k(a) = x^{k-1} a(1/x)$ for the reversal of a polynomial viewed with length $k$. If $A = Qf + R$ with $\deg R < d$ and $A$ of length $\ell$, then reversing gives

so the quotient is

a truncated product with the power series inverse of the reversal of $f$, and the remainder is $R = A - Qf$, of which only the low $d$ coefficients need to be computed. The reversal has constant term $\operatorname{lc}(f)$, which is why the leading coefficient must be invertible.

FLINT calls this inverse finv (or Binv, poly3inv, ...) and computes it with

nmod_poly_reverse(finv, f, f->length);
nmod_poly_inv_series(finv, finv, f->length);

A product of two reduced elements has length at most $2d - 1$, so the quotient has length at most $d$ and finv must be correct modulo $x^d$; computing it modulo $x^{d+1}$ as above is the convention. Every function with preinv in its name expects this finv: nmod_poly_mulmod_preinv, nmod_poly_divrem_newton_n_preinv, nmod_poly_div_newton_n_preinv, the powmod functions with preinv, nmod_poly_powers_mod_bsgs (which computes it itself) and the Brent–Kung modular composition functions. A reduction then costs two products of length about $d$, against the inversion and two products of nmod_poly_rem.

For moduli with at most five nonzero coefficients below the leading one, such as the trinomials and pentanomials returned by nmod_poly_minimal_irreducible, _nmod_poly_divrem_newton_n_preinv detects the sparsity (for $d > 20$) and reduces directly, with a few multiply-adds per coefficient of the dividend and no products.

Which function to use#

The preinv and precond functions require their inputs to be reduced modulo $f$; the plain mulmod and powmod functions reduce them first. nmod_poly_mulmod_precond precomputes either the matrix of multiplication by $a$ (fastest for $d$ up to a few hundred) or the quotient $\lfloor a x^d / f \rfloor$ (Shoup's method), and nmod_poly_mulmod_precond_init_num chooses between them from the expected number of multiplications. The precomputed object keeps pointers to $a$, $f$ and finv, which must not change while it is used.

The function $x \mapsto x^p$ is the Frobenius automorphism of $\mathbb{F}_{p^d}$, and $x^{p^i} \bmod f$ for successive $i$ is the core of distinct-degree factorisation. Once $h = x^p \bmod f$ is known, $x^{p^{i+1}} = (x^{p^i})^p = h(x^{p^i}) \bmod f$ is a modular composition, which Brent and Kung's algorithm computes in $O(d^{1/2} M(d) + d^{(\omega+1)/2})$ operations instead of the $O(M(d) \log p)$ of an exponentiation [GatSho1992] [KalSho1998]. The matrix of powers of $h$ can be precomputed with nmod_poly_precompute_matrix and reused for every step; see Modular composition.

Example: arithmetic in a finite field#

The program constructs $\mathbb{F}_{101^8}$ with the sparsest irreducible polynomial of degree 8 and checks some identities: the Frobenius automorphism has order 8, the multiplicative group has order $p^8 - 1$, inverses multiply to 1, and 1000 multiplications by a fixed element with a preconditioned object agree with one exponentiation.

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/nmod_poly.h>

int main(void)
{
    nmod_poly_t f, finv, a, b, t;
    nmod_poly_mulmod_precond_t pre;
    fmpz_t e;
    slong i;

    nmod_poly_init(f, 101);
    nmod_poly_init_mod(finv, f->mod);
    nmod_poly_init_mod(a, f->mod);
    nmod_poly_init_mod(b, f->mod);
    nmod_poly_init_mod(t, f->mod);
    fmpz_init(e);

    nmod_poly_minimal_irreducible(f, 8);
    flint_printf("f = "); nmod_poly_print_pretty(f, "x"); flint_printf("\n");

    /* finv = 1 / rev(f) mod x^(d+1) */
    nmod_poly_reverse(finv, f, f->length);
    nmod_poly_inv_series(finv, finv, f->length);

    /* x^(p^8) = x in GF(p^8), but x^(p^4) != x */
    fmpz_ui_pow_ui(e, 101, 8);
    nmod_poly_powmod_x_fmpz_preinv(t, e, f, finv);
    flint_printf("x^(p^8) == x: %d\n", nmod_poly_is_gen(t));
    fmpz_ui_pow_ui(e, 101, 4);
    nmod_poly_powmod_x_fmpz_preinv(t, e, f, finv);
    flint_printf("x^(p^4) == x: %d\n", nmod_poly_is_gen(t));

    /* a^(p^8 - 1) = 1 for a = x + 2 */
    nmod_poly_set_coeff_ui(a, 1, 1);
    nmod_poly_set_coeff_ui(a, 0, 2);
    fmpz_ui_pow_ui(e, 101, 8);
    fmpz_sub_ui(e, e, 1);
    nmod_poly_powmod_fmpz_binexp_preinv(t, a, e, f, finv);
    flint_printf("a^(p^8 - 1) == 1: %d\n", nmod_poly_is_one(t));

    /* inverse */
    nmod_poly_invmod(b, a, f);
    nmod_poly_mulmod_preinv(t, a, b, f, finv);
    flint_printf("a * a^-1 == 1: %d\n", nmod_poly_is_one(t));

    /* 1000 multiplications by a with a preconditioned object */
    nmod_poly_mulmod_precond_init_num(pre, a, f, finv, 1000);
    nmod_poly_one(b);
    for (i = 0; i < 1000; i++)
        nmod_poly_mulmod_precond(b, pre, b);
    nmod_poly_mulmod_precond_clear(pre);
    nmod_poly_powmod_ui_binexp_preinv(t, a, 1000, f, finv);
    flint_printf("a^1000 agrees: %d\n", nmod_poly_equal(b, t));

    fmpz_clear(e);
    nmod_poly_clear(f);
    nmod_poly_clear(finv);
    nmod_poly_clear(a);
    nmod_poly_clear(b);
    nmod_poly_clear(t);
    return 0;
}

f = x^8+2
x^(p^8) == x: 1
x^(p^4) == x: 0
a^(p^8 - 1) == 1: 1
a * a^-1 == 1: 1
a^1000 agrees: 1

The modulus is a binomial: an irreducible $x^8 + c$ exists because the only prime factor 2 of 8 divides $p - 1 = 100$ and, as $4 \mid 8$, also $4 \mid p - 1$. Reduction modulo $x^8 + 2$ costs one multiplication per coefficient, and products modulo $f$ are then dominated by the multiplication itself.