Composition and Taylor shifts

The composition $f(g)$ of polynomials of degrees $d$ and $e$ has degree $de$, so even writing it down costs $de$ operations. nmod_poly_compose computes it in $O(M(de) \log d)$ operations, close to this bound. The most common special case, $g = x + c$, is a Taylor shift, which changes the basis from powers of $x$ to powers of $x - c$ and costs only one multiplication of length $d$. For compositions modulo a third polynomial see Modular composition, and for power series composition nmod_poly_compose_series.

Algorithms#

nmod_poly_compose handles constant $g$ by an evaluation, and $g = a x^k + c$ by a Taylor shift followed by scaling the coefficients by powers of $a$ and spreading them out. Otherwise it uses Horner's rule, $f(g) = (\cdots(f_d g + f_{d-1}) g + \cdots) g + f_0$, for $d \le 6$, and the divide-and-conquer algorithm for larger $d$: with $f = f_{\mathrm{lo}} + x^{k} f_{\mathrm{hi}}$ for $k$ a power of two about $d/2$,

where the powers $g^{2^j}$ are computed once. The last level multiplies polynomials of length about $de/2$, which dominates the cost.

The Taylor shift $f(x + c)$ by Horner's rule takes $d^2/2$ multiply-adds (only additions for $c = \pm 1$) and works for every modulus. The convolution method [ASU1975] [GG1997] uses

in which the inner sums are the coefficients of one product of length $d + 1$, after reversing $(k! f_k)$. It divides by $d!$, so it needs a prime modulus larger than $d$ (or one without prime factors up to $d$); for other moduli it throws FLINT_IMPINV. nmod_poly_taylor_shift uses Horner's rule for lengths below 100, for lengths larger than the modulus, and for $c = \pm 1$ below length 1000, and the convolution method otherwise. A composite modulus with a small prime factor and a length between 100 and the modulus therefore needs nmod_poly_taylor_shift_horner.

Examples#

Composition agrees with evaluation, $f(g)(c) = f(g(c))$, and the Taylor shift by $c$ moves the roots by $-c$:

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

int main(void)
{
    flint_rand_t state;
    nmod_poly_t f, g, h;
    ulong p = 1000003, c = 12345, roots[3] = { 1, 2, 3 };

    flint_rand_init(state);
    nmod_poly_init(f, p);
    nmod_poly_init(g, p);
    nmod_poly_init(h, p);

    nmod_poly_rand(f, state, 200);
    nmod_poly_rand(g, state, 30);
    nmod_poly_compose(h, f, g);
    flint_printf("deg f(g) = %wd\n", nmod_poly_degree(h));
    flint_printf("f(g)(c) == f(g(c)): %d\n",
        nmod_poly_evaluate_nmod(h, c) ==
        nmod_poly_evaluate_nmod(f, nmod_poly_evaluate_nmod(g, c)));

    /* (x - 1)(x - 2)(x - 3) shifted by 10 has roots -9, -8, -7 */
    nmod_poly_product_roots_nmod_vec(f, roots, 3);
    nmod_poly_taylor_shift(g, f, 10);
    flint_printf("f(x + 10) = "); nmod_poly_print_pretty(g, "x");
    flint_printf("\ng(-9) = %wu\n", nmod_poly_evaluate_nmod(g, p - 9));

    /* the two Taylor shift algorithms agree */
    nmod_poly_rand(f, state, 5000);
    nmod_poly_taylor_shift_horner(g, f, c);
    nmod_poly_taylor_shift_convolution(h, f, c);
    flint_printf("horner == convolution: %d\n", nmod_poly_equal(g, h));

    nmod_poly_clear(f);
    nmod_poly_clear(g);
    nmod_poly_clear(h);
    flint_rand_clear(state);
    return 0;
}

deg f(g) = 5771
f(g)(c) == f(g(c)): 1
f(x + 10) = x^3+24*x^2+191*x+504
g(-9) = 0
horner == convolution: 1