Composition and Taylor shifts

fmpz_poly_compose computes $g(h(x))$, a polynomial of degree $\deg g \cdot \deg h$. The Taylor shift $f(x + c)$ is the special case $h = x + c$, which occurs often enough (for example in root isolation) to have its own functions. For composition of power series, see fmpz_poly_compose_series.

Algorithms#

Horner's rule $g(h) = (\cdots(g_m h + g_{m-1}) h + \cdots) h + g_0$ performs $m = \deg g$ products of increasingly unbalanced size. The divide-and-conquer algorithm groups the coefficients of $g$ in pairs, $g_{2i} + g_{2i+1} h$, then combines pairs with $h^2$, quadruples with $h^4$ and so on, so that all products are balanced; with fast multiplication its cost is $O(M(mn) \log m)$ for $n = \deg h$, against $O(m M(mn))$ for Horner's rule. fmpz_poly_compose uses Horner's rule for $\operatorname{len}(g) \le 4$ and divide and conquer otherwise, and reduces inner polynomials of the form $a x^k + b$ to a Taylor shift, a scaling and an inflation.

For the Taylor shift, Horner's scheme $f \mapsto f(x + c)$ costs $O(n^2)$ additions and multiplications by $c$ (only additions for $c = \pm 1$), on coefficients that grow to $O(n \log |c| + n + b)$ bits. The divide-and-conquer algorithm writes $f = f_0 + x^m f_1$, shifts the halves recursively and computes $(x + c)^m$ by the binomial theorem, so that the work is in a few large products [GG1997]. fmpz_poly_taylor_shift uses Horner's scheme below length 64 and divide and conquer otherwise; the divide-and-conquer function itself returns to Horner's scheme below a length between 100 and 1000 that grows with the coefficient size, and runs the two recursive halves in parallel when threads are available. fmpz_poly_taylor_shift_multi_mod computes the shift modulo many word-size primes in parallel and reconstructs the result, which is competitive for long polynomials with small coefficients and many threads.

Examples#

Composition with $h = x^2 + 1$ and a Taylor shift by $-1$; the shift of a polynomial by $c$ and then by $-c$ returns the original:

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

int main(void)
{
    fmpz_poly_t g, h, r, s;
    fmpz_t c;

    fmpz_poly_init(g);
    fmpz_poly_init(h);
    fmpz_poly_init(r);
    fmpz_poly_init(s);
    fmpz_init(c);

    fmpz_poly_set_str(g, "3  1 -3 2");        /* 2x^2 - 3x + 1 */
    fmpz_poly_set_str(h, "3  1 0 1");         /* x^2 + 1 */
    fmpz_poly_compose(r, g, h);
    flint_printf("g(h(x)) = "); fmpz_poly_print_pretty(r, "x"); flint_printf("\n");

    fmpz_set_si(c, -1);
    fmpz_poly_taylor_shift(r, g, c);
    flint_printf("g(x - 1) = "); fmpz_poly_print_pretty(r, "x"); flint_printf("\n");

    /* the shift agrees with composition by x + c */
    fmpz_poly_chebyshev_t(g, 100);
    fmpz_set_si(c, 7);
    fmpz_poly_taylor_shift(r, g, c);
    fmpz_poly_set_str(h, "2  7 1");
    fmpz_poly_compose(s, g, h);
    flint_printf("shift == compose: %d\n", fmpz_poly_equal(r, s));
    fmpz_neg(c, c);
    fmpz_poly_taylor_shift(r, r, c);
    flint_printf("shift back: %d\n", fmpz_poly_equal(r, g));

    fmpz_poly_clear(g);
    fmpz_poly_clear(h);
    fmpz_poly_clear(r);
    fmpz_poly_clear(s);
    fmpz_clear(c);
    return 0;
}

g(h(x)) = 2*x^4+x^2
g(x - 1) = 2*x^2-7*x+6
shift == compose: 1
shift back: 1