Composition and modular composition

Composition f(g) of polynomials appears in three forms: exact composition of polynomials, composition of power series truncated to x^n (which requires g(0) = 0), and modular composition f(g) \bmod h, the basic operation of polynomial factoring over finite fields [GatSho1992; KalSho1998]. In all three, Horner's rule f(g) = f_0 + g(f_1 + g(f_2 + \cdots)) uses \ell_1 - 1 multiplications by g, which is too slow for long f. The faster methods split f into blocks and reuse powers of g, or (for power series) reorganise the computation entirely. Brent and Kung [BrentKung1978] introduced the baby-step giant-step method, which is also the basis of evaluating polynomials at matrices by Paterson and Stockmeyer [PS1973]. Kinoshita and Li [KL2024] recently found a quasi-linear algorithm for power series composition.

Throughout, f has length \ell_1, g has length \ell_2, M(n) is the cost of multiplying two polynomials of length n, and \omega is an exponent of matrix multiplication (\omega \le \log_2 7 by Strassen's method).

Generic version#

All the composition algorithms described here are implemented once in gr_poly for any commutative ring, and nmod_poly, fmpz_poly, fmpz_mod_poly, fq_poly, arb_poly and acb_poly call them with their contexts. The generic dispatchers choose by size only, as described below for nmod_poly (over balls, the cutoff between Brent–Kung and Kinoshita–Li depends on the precision). The Taylor shift, a ring method, uses Horner's scheme up to length 20 and divide and conquer above, with fmpz and the ball rings supplying their own versions. Modular composition with the Brent-Kung algorithm multiplies matrices with gr_mat_mul, so it uses the fastest matrix multiplication of the ring (Linear algebra over generic rings and its specialisations). Series reversion by Lagrange inversion divides by integers and fails in small characteristic, where Newton iteration is used. The same algorithms over other rings are summarised in Composition and Taylor shifts.

Composition of polynomials#

The result f(g) has length (\ell_1 - 1)(\ell_2 - 1) + 1. Horner's rule multiplies polynomials of growing length i(\ell_2 - 1) + 1 by g, which costs about as much as \tfrac{1}{2}\ell_1^2 products of length \ell_2. Divide and conquer instead pairs up the coefficients:

which is a polynomial of half the length in g^2, with polynomial "coefficients" f_{2i} + f_{2i+1} g. Repeating with g^2, g^4, \ldots gives:

Taylor shifts#

The special case g = x + c, the Taylor shift f(x + c), has a quasi-linear algorithm [ASU1975]; von zur Gathen and Gerhard [GG1997] compare this and other methods.

More generally, g = a x^d + c is a Taylor shift by c, a scaling of the coefficients by powers of a and the substitution x \to x^d.

In FLINT. nmod_poly_compose uses Horner's rule (nmod_poly_compose_horner) for \ell_1 \le 7. Otherwise it calls _gr_poly_compose, which handles g = a x^d + c by a Taylor shift, scaling and inflation, uses Horner's rule for \ell_1 \le 7, and Algorithm 1 otherwise. nmod_poly_taylor_shift uses Horner's rule (O(\ell_1^2) additions and multiplications) for lengths below 100, for lengths exceeding the modulus and for c = \pm 1 below length 1000, and Lemma 2 (nmod_poly_taylor_shift_convolution) otherwise. See also Composition and Taylor shifts.

Composition of power series#

Let g(0) = 0 and n the target precision. Horner's rule and divide and conquer work as before, truncating every product to length n. Both are fast when the exact composition is not much longer than n, but in general they need \Theta(\ell_1) products of length n.

Brent–Kung#

Let m \approx \sqrt{n} and split f into blocks of m coefficients:

The baby steps are g^0, g^1, \ldots, g^{m-1} \bmod x^n. Then all F_i(g) together are the rows of a matrix product: the matrix of coefficients of the F_i (about \ell_1/m rows and m columns) times the matrix whose rows are the baby steps (m rows, n columns). The giant steps combine the F_i(g) by Horner's rule in g^m.

Kinoshita–Li#

The algorithm of Kinoshita and Li [KL2024] rests on the identity, valid for \ell_1 \le n,

which follows by expanding 1/(1 - yg) = \sum_j y^j g^j. Writing Q(x, y) = 1 - y g(x), the bivariate rational function is simplified by Graeffe steps: Q(x, y) Q(-x, y) = V(x^2, y), so

which halves the x-degree that matters while the y-degree of the denominator doubles. After O(\log n) steps the denominator no longer depends on x; the problem is then a univariate power series division, and the result is lifted back through the stored Q's, keeping at each level only the slice of y-coefficients that contributes to [y^{n-1}]. Each level costs O(M(n)) with bivariate products done by Kronecker substitution, for O(M(n) \log n) in total.

In FLINT. nmod_poly_compose_series calls _gr_poly_compose_series, which uses direct formulas for n \le 3 and for monomials g = c x^d, Horner's rule if \ell_1 \le 7 or n \le 7, divide and conquer if \ell_1^2 < n or (\ell_1 - 1)(\ell_2 - 1) + 1 < 4n (when the exact composition is not much longer than n), Brent–Kung for n < 500 and Kinoshita–Li (_gr_poly_compose_series_kinoshita_li) for n \ge 500.

Reversion. The compositional inverse of f (with f(0) = 0 and f'(0) a unit), nmod_poly_revert_series, uses Johansson's fast Lagrange inversion [Joh2015b], a baby-step giant-step method, for n < 2500, and Newton iteration on f(g) = x, which needs one composition and one series division per step, for larger n or when Lagrange inversion fails because some k < n is not invertible.

Modular composition#

Let h have degree n with invertible leading coefficient and let f, g have degree less than n. Horner's rule computes f(g) \bmod h with \ell_1 - 1 multiplications modulo h, each costing O(M(n)) with Newton division and a precomputed inverse of \operatorname{rev}(h) (Division with remainder): O(n M(n)) in total. The Brent–Kung method works modulo h just as modulo x^n:

Correctness follows from f(g) = \sum_i F_i(g) (g^m)^i and the fact that reduction modulo h commutes with sums and products. Since m^2 > n \ge \ell_1, the matrix B holds all coefficients of f.

With classical matrix multiplication the matrix step costs n^2 operations, as much as schoolbook arithmetic; the method pays off because matrix multiplication over \mathbb{Z}/p\mathbb{Z} is very efficient in practice (Strassen, BLAS and delayed reduction, see Multiplying matrices over Z/nZ), and because the matrix A depends only on g and h, so it can be reused for many f. Nüsken and Ziegler [NZ2004] improve the exponent using rectangular matrix multiplication, and Kedlaya and Umans [KU2011] give an algorithm of nearly linear complexity over finite fields; neither is implemented in FLINT.

FLINT improves step 4 with rectangular splitting when the matrix has more than 10 rows: it precomputes H^0, \ldots, H^t for t \approx \sqrt{m} (at most 30), accumulates each block \sum_{j < t} C_{it + j} H^j without reduction (each term is a product of two polynomials of degree less than n), reduces once per block, and combines the blocks by Horner's rule in H^t. This replaces m modular multiplications by about m plain products and m/t + t reductions.

In FLINT.

• nmod_poly_compose_mod reduces g modulo h and uses Horner's rule (nmod_poly_compose_mod_horner) if \deg h < 7 or \ell_1 \ge n + 1, and Brent–Kung (nmod_poly_compose_mod_brent_kung, which requires \ell_1 \le n) otherwise. The latter calls _gr_poly_compose_mod_brent_kung, which computes the baby steps as g^i = g^{\lceil i/2 \rceil} g^{\lfloor i/2 \rfloor} and multiplies the matrices with gr_mat_mul.
• nmod_poly_compose_mod_brent_kung_preinv takes the inverse of the reversal of h as an extra argument; it computes the baby steps with _nmod_poly_powers_mod_preinv_naive, the matrix product with nmod_mat_mul, H = g^{\lfloor m/2 \rfloor} g^{m - \lfloor m/2 \rfloor} and the giant steps with the rectangular splitting described above (_nmod_poly_mod_matrix_rows_evaluate).
• nmod_poly_precompute_matrix computes the matrix A of baby steps once, and nmod_poly_compose_mod_brent_kung_precomp_preinv reuses it for several f. The distinct-degree and equal-degree factorisation and irreducibility tests of gr_poly use the same technique to compose many polynomials with the same g modulo the same h [KalSho1998].
• nmod_poly_powers_mod_bsgs computes f^0, \ldots, f^{k-1} \bmod g by a baby-step giant-step decomposition (\sqrt{k} baby steps, then giant steps, then the remaining powers in parallel blocks on the thread pool); nmod_poly_powers_mod_naive computes them one after the other.

Modular composition describes these functions, with Rabin's irreducibility test as an example.