vii. Composition and Taylor Shift

$f(g(x))\mathrm{mod}{x}^n$, the compositional inverse (reversion) $h$ with $g(h(x))=x$, power projection ${\sum }_i{y}^i[x^{n-1}]w{g}^i$, Lagrange inversion, and the Taylor shift $f(x+c)$ with its applications (sample shifting, ${\sum }_i\left(\genfrac{}{}{0ex}{}{a_i}{k}\right)$ for all $k$). Uses the core of NTT and Power Series, and shiftSamples from Evaluation and Interpolation.

Operation	Method	Cost
$f(g)\mathrm{mod}{x}^n$	Brent–Kung baby-step giant-step	$O(n\sqrt{n}\log n+n^2)$, fine for $n\le 2\cdot 10^4$
$f(g)\mathrm{mod}{x}^n$	Kinoshita–Li (transposed power projection)	$O(n{\log }^{2}n)$ ($n=10^5$ in about 2 s)
power projection $a_i=[x^{n-1}]w{g}^i$, $i<m$	Graeffe / Bostan–Mori in two variables on $\frac{w(x)}{1-yg(x)}$	$O(n{\log }^{2}n)$
reversion $g(h(x))=x$	Newton $h\leftarrow h-\frac{g(h)-x}{g^{\prime }(h)}$ with composition	$O(n^2)$ with Brent–Kung
reversion	power projection + Lagrange inversion + $N$-th root	$O(n{\log }^{2}n)$
single coefficient $[x^N]h^k$ of an inverse	Lagrange: $N[x^N]h^k=k[x^{N-k}](x/g{)}^N$	one power, $O(n\log n)$
$f(x+c)$	one convolution	$O(n\log n)$
$f(m),\dots ,f(m+cnt-1)$ from $f(\mathrm{0..}n-1)$	sample shift (convolution)	$O((n+cnt)\log )$
${\sum }_i\left(\genfrac{}{}{0ex}{}{a_i}{k}\right)$ for all $k$	Taylor shift by 1 of ${\sum }_i{x}^{a_i}$	$O(A\log A)$

Poly compose(const Poly& f, const Poly& g, int n) {        // f(g(x)) mod x^n, any g (Brent-Kung)
    Poly res(n, 0); if (f.empty() || n == 0) return res;
    int d = sz(f), k = max(1, (int)sqrt((double)d)); Poly gm = polyTrunc(g, n);
    vector<Poly> pw(k + 1); pw[0] = Poly(n, 0); pw[0][0] = 1;          // g^0 .. g^k mod x^n
    for (int i = 1; i <= k; i++) pw[i] = mulTrunc(pw[i - 1], gm, n);
    int blocks = (d + k - 1) / k;
    for (int j = blocks - 1; j >= 0; j--) {                 // Horner in g^k: res = res g^k + F_j(g)
        if (j != blocks - 1) res = mulTrunc(res, pw[k], n);
        for (int i = 0; i < k && j * k + i < d; i++) { mint c = f[j * k + i]; if (c == mint(0)) continue; for (int t = 0; t < n; t++) res[t] += c * pw[i][t]; }
    }
    return res;
}
Poly compositionalInverse(const Poly& g, int n) {           // g(0) = 0, g'(0) != 0; Newton, O(n^2)
    Poly h = {0, g[1].inv()}, dg = derivative(g);
    for (int len = 2; len < n;) {
        len = min(2 * len, n);
        Poly gh = compose(g, h, len); gh[1] -= 1;           // g(h) - x
        Poly step = mulTrunc(gh, polyInv(compose(dg, h, len), len), len);
        h.resize(len); for (int i = 0; i < len; i++) h[i] -= step[i];
    }
    h.resize(n); return h;
}
Poly taylorShift(const Poly& f, mint c) {                   // f(x + c), O(n log n)
    int n = sz(f); if (n == 0) return f;
    vector<mint> fact(n, 1), ifact(n); for (int i = 1; i < n; i++) fact[i] = fact[i - 1] * mint(i);
    ifact[n - 1] = fact[n - 1].inv(); for (int i = n - 1; i > 0; i--) ifact[i - 1] = ifact[i] * mint(i);
    Poly a(n), b(n); mint cp = 1;
    for (int i = 0; i < n; i++) a[n - 1 - i] = f[i] * fact[i], b[i] = cp * ifact[i], cp *= c;
    Poly prod = polyMul(a, b), g(n);
    for (int k = 0; k < n; k++) g[k] = prod[n - 1 - k] * ifact[k];
    return g;
}
vector<mint> samplingShift(const vector<mint>& ys, ll m, int cnt, const Comb& cb) {   // ys[i] = f(i), deg f < n; f(m..m+cnt-1)
    int n = sz(ys); vector<mint> res; if (cnt <= 0) return res;
    if (n == 0) return vector<mint>(cnt, 0);
    while (m < n && sz(res) < cnt) res.push_back(ys[m++]);   // known values
    if (sz(res) < cnt) { auto rest = shiftSamples(ys, mint(m), cnt - sz(res), cb); res.insert(res.end(), all(rest)); }
    return res;
}
vector<mint> binomialSumsAll(const vector<int>& a, int K) { // s[k] = sum_i C(a_i, k), a_i >= 0
    int mx = 0; for (int x : a) mx = max(mx, x);
    Poly p(mx + 1, 0); for (int x : a) p[x] += 1;
    Poly q = taylorShift(p, 1); vector<mint> res(K + 1, 0);
    for (int k = 0; k <= K && k <= mx; k++) res[k] = q[k];
    return res;
}

Proof: Brent–Kung, Newton reversion, Taylor shift, binomial sums

Brent–Kung. Write $f={\sum }_j{F}_j(x)x^{jk}$ with $\mathrm{deg}{F}_j<k$; then $f(g)={\sum }_j{F}_j(g)(g^k{)}^j$, and each $F_j(g)$ is a linear combination of the stored $g^0..g^{k-1}$ ($O(kn)$ per block, $O(dn)$ total), combined by Horner in $g^k$ with $d/k$ truncated multiplications: $O(dn+\frac{d}{k}n\log n)$, minimised near $k=\sqrt{d}$. Reversion: solve $G(h)=g(h)-x=0$; Newton’s step for a composition is $h-\frac{g(h)-x}{g^{\prime }(h)}$ ($\frac{d}{dh}g(h)=g^{\prime }(h)$), doubling the correct coefficients from $h=g_1^{-1}x$ (correct mod $x^2$). Taylor shift: $[x^k]f(x+c)={\sum }_{i\ge k}{f}_i\left(\genfrac{}{}{0ex}{}{i}{k}\right)c^{i-k}=\frac{1}{k!}{\sum }_i(f_{i}i!)\frac{c^{i-k}}{(i-k)!}$, a convolution of the reversed $(f_{i}i!)$ with $(c^j/j!)$. Binomial sums: ${\sum }_i(1+x{)}^{a_i}$ is ${\sum }_i{x}^{a_i}$ shifted by 1, and $[x^k](1+x{)}^a=(\genfrac{}{}{0ex}{}{a}{k})$. Sample shift is proved in Evaluation and Interpolation. All checked against naive composition and direct evaluation.

Kinoshita–Li: power projection and fast composition

Power projection. ${\sum }_{i<m}{y}^i[x^{n-1}]w(x)g(x{)}^i=[x^{n-1}]\frac{w(x)}{1-yg(x)}$ as a series in $y$ truncated at $y^m$. Bostan–Mori in two variables: multiply numerator and denominator by $Q(-x,y)$; the denominator $Q(x,y)Q(-x,y)=V(x^2,y)$ is even in $x$, so $[x^N]\frac{P(x,y)}{V(x^2,y)}=[x^{\lfloor N/2\rfloor }]\frac{P_{N\mathrm{mod}2}(x,y)}{V(x,y)}$ where $P_r$ keeps the coefficients of $x^{2k+r}$. Each round halves the $x$-degree and doubles the $y$-degree (truncated at $m$), so all rounds cost $O(n{\log }^{2}n)$; at $N=0$ the answer is $\frac{P(0,y)}{Q(0,y)}\mathrm{mod}{y}^m$. Composition is the transpose: $⟨w,f(g)\mathrm{mod}{x}^n⟩$ (pairing $w$ reversed with $f(g)$) equals ${\sum }_i{f}_i[x^{n-1}]w{g}^i$, so the linear map $f\mapsto f(g)\mathrm{mod}{x}^n$ is the transpose of $w\mapsto$ (power projection). Running the transposed steps of power projection in reverse order (Tellegen’s principle) computes composition in the same $O(n{\log }^{2}n)$. Reversion: Lagrange inversion turns $p_k=[x^N]g^k$ (one power projection) into $(x/h{)}^N\mathrm{mod}{x}^N$, then an $N$-th root by $\log$ / $\exp$.

vector<mint> mul2D(const vector<mint>& a, int ax, int ay, const vector<mint>& b, int bx, int by) {   // row-major x * Y + y
    int cy = ay + by - 1, cx = ax + bx - 1;
    Poly fa((ax - 1) * cy + ay, 0), fb((bx - 1) * cy + by, 0);
    for (int x = 0; x < ax; x++) for (int y = 0; y < ay; y++) fa[x * cy + y] = a[x * ay + y];
    for (int x = 0; x < bx; x++) for (int y = 0; y < by; y++) fb[x * cy + y] = b[x * by + y];
    Poly c = polyMul(fa, fb); c.resize(cx * cy, 0); return c;
}
struct GraeffeLevels {                                      // Q_0 = 1 - y g(x), Q_{t+1}(x^2, y) = Q_t(x, y) Q_t(-x, y), truncated
    vector<vector<mint>> Q; vector<int> rows, cols, N;
    GraeffeLevels(const Poly& g, int n, int m) {
        int qy = min(2, m), qx = n, deg = n - 1; vector<mint> q(qx * qy, 0); q[0] = 1;
        if (qy == 2) for (int x = 0; x < qx; x++) q[x * 2 + 1] = -(x < sz(g) ? g[x] : mint(0));
        while (true) {
            Q.push_back(q), rows.push_back(qx), cols.push_back(qy), N.push_back(deg);
            if (deg == 0) break;
            vector<mint> qn = q; for (int x = 1; x < qx; x += 2) for (int y = 0; y < qy; y++) qn[x * qy + y] = -qn[x * qy + y];
            auto prod = mul2D(q, qx, qy, qn, qx, qy);
            int py = 2 * qy - 1, nd = deg / 2, nqx = nd + 1, nqy = min(py, m); vector<mint> v(nqx * nqy);
            for (int k = 0; k < nqx; k++) for (int y = 0; y < nqy; y++) v[k * nqy + y] = prod[(2 * k) * py + y];
            q = v, qx = nqx, qy = nqy, deg = nd;
        }
    }
};
vector<mint> powerProjection(const Poly& w, const Poly& g, int m) {   // a_i = [x^(|w|-1)] w g^i, i < m
    int n = sz(w); if (n == 0 || m == 0) return vector<mint>(m, 0);
    GraeffeLevels lv(g, n, m); vector<mint> p(all(w)); int px = n, py = 1;
    for (int t = 0; t + 1 < sz(lv.Q); t++) {
        int qx = lv.rows[t], qy = lv.cols[t], deg = lv.N[t];
        vector<mint> qn = lv.Q[t]; for (int x = 1; x < qx; x += 2) for (int y = 0; y < qy; y++) qn[x * qy + y] = -qn[x * qy + y];
        auto u = mul2D(p, px, py, qn, qx, qy);
        int uy = py + qy - 1, par = deg & 1, nx = deg / 2 + 1, ny = min(uy, m); vector<mint> np(nx * ny);
        for (int k = 0; k < nx; k++) for (int y = 0; y < ny; y++) np[k * ny + y] = u[(2 * k + par) * uy + y];
        p = np, px = nx, py = ny;
    }
    int qy = lv.cols.back();
    Poly q0(lv.Q.back().begin(), lv.Q.back().begin() + qy), p0(p.begin(), p.begin() + py);
    return mulTrunc(p0, polyInv(q0, m), m);
}
Poly composeFast(Poly f, Poly g, int n) {                   // f(g(x)) mod x^n, O(n log^2 n)
    if (n == 0) return {};
    g = polyTrunc(g, n);
    if (g[0] != mint(0)) f = taylorShift(f, g[0]), g[0] = 0;
    if (sz(f) > n) f.resize(n);                             // g(0) = 0: g^i = O(x^i)
    int m = sz(f); if (m == 0) return Poly(n, 0);
    GraeffeLevels lv(g, n, m); int T = sz(lv.Q) - 1;
    vector<int> D(T + 1, 0); for (int t = 0; t < T; t++) D[t + 1] = min(D[t] + lv.cols[t] - 1, m - 1);   // y-degree bounds
    int qy = lv.cols.back(); Poly q0(lv.Q.back().begin(), lv.Q.back().begin() + qy);
    Poly R = polyInv(q0, m), Rr(R.rbegin(), R.rend()), fb = polyMul(f, Rr);
    vector<mint> B(D[T] + 1); for (int j = 0; j <= D[T]; j++) B[j] = j + m - 1 < sz(fb) ? fb[j + m - 1] : mint(0);
    int by = D[T] + 1;
    for (int t = T - 1; t >= 0; t--) {                       // transposed rounds, last to first
        int qx = lv.rows[t], qyt = lv.cols[t], deg = lv.N[t], par = deg & 1, ey = D[t + 1] + 1, nx = deg / 2 + 1;
        vector<mint> E(qx * ey, 0);
        for (int k = 0; k < nx; k++) for (int y = 0; y < ey && y < by; y++) E[(2 * k + par) * ey + y] = B[k * by + y];
        vector<mint> qr(qx * qyt);                          // Q_t(-x, y) reversed in both variables
        for (int c = 0; c < qx; c++) for (int d = 0; d < qyt; d++) { mint v = lv.Q[t][c * qyt + d]; if (c & 1) v = -v; qr[(qx - 1 - c) * qyt + (qyt - 1 - d)] = v; }
        auto conv = mul2D(E, qx, ey, qr, qx, qyt);
        int cy = ey + qyt - 1, nb = D[t] + 1; vector<mint> nbv(qx * nb);
        for (int a = 0; a < qx; a++) for (int b = 0; b < nb; b++) nbv[a * nb + b] = conv[(a + qx - 1) * cy + (b + qyt - 1)];
        B = nbv, by = nb;
    }
    Poly h(n); for (int j = 0; j < n; j++) h[j] = B[(n - 1 - j) * by];
    return h;
}
Poly compositionalInverseFast(const Poly& g, int n) {       // g(0) = 0, g'(0) != 0, O(n log^2 n)
    if (n <= 1) return Poly(n, 0);
    int N = n - 1; Poly w(n, 0); w[0] = 1;
    auto p = powerProjection(w, polyTrunc(g, n), N + 1);    // p_k = [x^N] g^k
    Poly S(N); for (int j = 0; j < N; j++) S[j] = mint(N) * mint(N - j).inv() * p[N - j];   // (x/h)^N mod x^N
    mint c = S[0].inv(); for (auto& x : S) x *= c;
    Poly lg = polyLog(S, N); mint iN = mint(N).inv(); for (auto& x : lg) x *= iN;
    Poly xh = polyExp(lg, N); for (auto& x : xh) x *= g[1];  // x / h, constant term g'(0)
    Poly inv = polyInv(xh, N), h(n, 0); for (int j = 0; j < N; j++) h[j + 1] = inv[j];
    return h;
}

Proof: Graeffe power projection, transposition and Lagrange inversion

Power projection. In the ring of polynomials in $x$ over formal power series in $y$, the quotient $\frac{P(x,y)}{Q(x,y)}$ equals $\frac{P(x,y)Q(-x,y)}{V(x^2,y)}$; since the denominator has only even powers of $x$, the coefficient of $x^N$ comes only from numerator terms $x^{2k+(N\mathrm{mod}2)}$, giving $[x^{\lfloor N/2\rfloor }]\frac{P_{N\mathrm{mod}2}(x,y)}{V(x,y)}$. Coefficients of $x$ above the current $N$ never matter and $y$-powers $\ge m$ are discarded (the answer is wanted mod $y^m$), so the sizes stay $O(n)$ per round as $N$ halves while the $y$-degree doubles, and each round is one 2D product: $O(n\log n)$ per round, $\log n$ rounds. When $N=0$, $[x^0]\frac{P}{Q}=\frac{P(0,y)}{Q(0,y)}$. With $Q_0=1-yg(x)$ and $P=w$, $\frac{w}{1-yg}={\sum }_i{y}^{i}w{g}^i$, which is the claim.

Composition as transpose. Both maps are linear: $C:f\mapsto f(g)\mathrm{mod}{x}^n$ and $\Pi :w\mapsto (a_i{)}_{i<m}$ with $a_i={\sum }_j{w}_j[x^{n-1-j}]g^i$. Then ${\sum }_i{v}_i{a}_i={\sum }_j{w}_j[x^{n-1-j}]{\sum }_i{v}_i{g}^i={\sum }_j{w}_{j}C(v{)}_{n-1-j}$, so $\Pi =J{C}^T$ with the index reversal $J$. Power projection is a composition of linear maps in $w$ (2D multiplications by fixed $Q_t(-x,y)$, extraction of parity classes and truncations, the final multiplication by $Q(0,y{)}^{-1}$); composeFast applies their transposes in reverse order: a transposed multiplication by a fixed polynomial is a multiplication by its reversal followed by extracting the middle coefficients (conv[(a + qx - 1)\cdot cy + (b + qyt - 1)]), the transpose of a parity extraction is the corresponding injection (E), and the first step is the transposed series multiplication f * Rr. Transposition preserves the cost (Tellegen’s principle). A non-zero $g(0)$ is first moved into $f$ by a Taylor shift. Verified against Brent–Kung on random inputs of many sizes.

Lagrange inversion. For $g(0)=0\ne g^{\prime }(0)$ with compositional inverse $h$ and $N,k\ge 1$: $N[x^N]g^k=k[x^{N-k}](x/h{)}^N$. Proof with formal residues ($\mathrm{res}A=[x^{-1}]A$ for Laurent series): $\mathrm{res}{A}^{\prime }=0$, and for $u$ of valuation 1, ${\mathrm{res}}_{x}A(x)={\mathrm{res}}_{t}A(u(t))u^{\prime }(t)$ (check on $x^j$: for $j\ne -1$ the right side is the residue of a derivative, for $j=-1$ it is $\mathrm{res}\frac{u^{\prime }}{u}=1$). Then $N[x^N]g^k={\mathrm{res}}_x{x}^{-N}(g^k{)}^{\prime }(x)$; substituting $x=h(t)$ gives ${\mathrm{res}}_t{h}^{-N}\frac{d}{dt}g(h(t){)}^k={\mathrm{res}}_t{h}^{-N}k{t}^{k-1}=k[t^{N-k}](t/h{)}^N$. So the projection $p_k=[x^N]g^k$ yields $[x^{N-k}](x/h{)}^N=\frac{N}{k}{p}_k$ for $k=\mathrm{1..}N$, i.e. $(x/h{)}^N\mathrm{mod}{x}^N$. Its constant term is $g^{\prime }(0{)}^N$; taking $\exp (\frac{1}{N}\log (\cdot /S_0))$ and multiplying by $g^{\prime }(0)$ (the constant of $x/h$) gives $x/h$ (the $N$-th root with that constant is unique), and $h=x\cdot (x/h{)}^{-1}$. Needs $N<p$. Verified: $g(h)=x$ for random $g$, against the Newton version.

Related notes

• Linear Recurrences (Bostan–Mori in one variable)
• Generating Functions (Lagrange inversion for trees)
• Evaluation and Interpolation