形式幂级数复合|复合逆

形式幂级数的复合和复合逆也是常见的形式幂级数操作，对于没有特殊性质的 $f$ 之前我们一直使用的多是 $O (n^{2})$ 的算法（该算法仍需使用 FFT）来计算 $f (g) mod x^{n}$ 其中 $f \in C [[x]], g \in x C [[x]]$ ，但是因为效率较低应用较少．我们介绍 Kinoshita–Li 的 $O (M (n) \log n)$ 的算法，其中 $O (M (n))$ 为两个次数为 $O (n)$ 的多项式相乘的时间．

形式幂级数/多项式的复合

若要计算 $f (g (x)) mod x^{n}$ 那么需要 $f (g (x))$ 的每一项系数都是有限项之和，所以之前要求 $f (x) \in C [[x]], g (x) \in x C [[x]]$ ，而如果 $f (x), g (x) \in C [x]$ 也可以满足这个条件．因为我们需要将 $f (g (x))$ 的系数截断，不妨直接考虑 $f (x), g (x)$ 都是多项式的情况．对于 $f (x) = \sum_{j = 0}^{n - 1} f_{j} x^{j}$ ，有

f (g (x)) = \sum_{j = 0}^{n - 1} f_{j} g {(x)}^{j}

我们考虑环 $C [x] ((y))$ 上的有理函数

\begin{aligned} \frac{f (y^{- 1})}{1 - y \cdot g (x)} & = \sum_{j \geq 0} (\dots + f_{j} y^{- j} + \dots) g (x)^{j} y^{j} \\ f (g (x)) & = [y^{0}] \frac{f (y^{- 1})}{1 - y \cdot g (x)} \end{aligned}

根据常系数齐次线性递推中提到的 Bostan–Mori 算法，Kinoshita 和 Li 指出可以将其修改为二元形式：

\begin{aligned} \frac{P (y)}{Q (x, y)} mod x^{n} & = (\frac{P (y)}{Q (x, y) Q (- x, y)} mod x^{n}) Q (- x, y) mod x^{n} \\ = (\frac{P (y)}{V (x^{2}, y)} mod x^{n}) Q (- x, y) mod x^{n} \\ = {(\frac{P (y)}{V (z, y)} mod z^{⌈ n / 2 ⌉}) |}_{z = x^{2}} Q (- x, y) mod x^{n} \end{aligned}

这样递归的计算在 $n = 1$ 时我们只需计算

\frac{P (y)}{Q (x, y)} mod x = \frac{P (y)}{Q (0, y)} \in C ((y))

在计算 $\frac{P (y)}{V (z, y)} mod z^{⌈ n / 2 ⌉} \in C [z] ((y))$ 时我们不需要保留所有 $y$ 的系数，因为最后我们只需要提取 $y^{0}$ 的系数，所以 $y^{> 0}$ 的系数是不需要的，而因为求出前者之后需要将其乘以若干个形如 $Q (- x, y) \in C [x, y]$ 的「多项式」，所以只需要保留对于 $y^{0}$ 有贡献的系数即可．我们准备好给出伪代码：

\begin{array}{ll} Algorithm Comp (P (y), Q (x, y), n, m) : \\ Input : P = \sum_{0 \leq j < n} p_{j} y^{- j} \in C ((y)), Q \in C [x, y], n, m \in N_{> 0} . \\ Output : [y^{(- m, 0]}] \frac{P (y)}{Q (x, y)} mod x^{n} . \\ Require : [x^{0} y^{0}] Q = 1 . \\ 1 & if n = 1 then return ([y^{- m + 1}] \frac{P (y)}{Q (0, y)}, \dots, [y^{0}] \frac{P (y)}{Q (0, y)}) \\ 2 & V (x^{2}, y) \leftarrow Q (x, y) Q (- x, y) mod x^{n} \\ 3 & d \leftarrow \deg_{y} Q (- x, y) \\ 4 & (t_{- (m + d) + 1}, \dots, t_{0}) \leftarrow Comp (P (y), V (x, y), ⌈ n / 2 ⌉, m + d) \\ 5 & T (x, y) \leftarrow \sum_{j = - (m + d) + 1}^{0} t_{j} y^{j} \\ 6 & U (x, y) = \sum_{j = - (m + d) + 1}^{d} u_{j} y^{j} \leftarrow T (x^{2}, y) Q (- x, y) mod x^{n} \\ 7 & return (u_{- m + 1}, \dots, u_{0}) \end{array}

那么我们有

f (g (x)) mod x^{n} = Comp (f (y^{- 1}), 1 - y \cdot g (x), max {1 + \deg f, n}, 1) mod x^{n}

注意第三个参数是因为 $g (0)$ 可能不为零，如果 $\deg f \geq n$ 此时不能截断 $f (x)$ 来计算 $f (g (x))$ ，我们也可以选择计算 $f (g) = f \circ (x + g (0)) \circ (g - g (0))$ ，此时可以取 $F := f (x + g (0)) mod x^{n}$ 和 $G := g - g (0)$ 转而计算 $Comp (F (y^{- 1}), 1 - y \cdot G (x), n, 1)$ ．

另外因为调用的限制最后递归终止时的 $Q (0, y)^{- 1}$ 是可以直接导出的，不需要使用形式幂级数的乘法逆元算法来计算，我们只需计算一次乘法然后提取需要的系数．

常见的特殊形式复合

我们常用的多项式初等函数都可以通过复合计算：

\begin{aligned} g (0) = 1 & , g^{- 1} = 1 + (1 - g) + (1 - g)^{2} + \dots \\ g (0) = 1 & , \log g = - \frac{1 - g}{1} - \frac{(1 - g)^{2}}{2} - \frac{(1 - g)^{3}}{3} - \dots \\ g (0) = 0 & , \exp g = 1 + \frac{g}{1!} + \frac{g^{2}}{2!} + \frac{g^{3}}{3!} + \dots \\ g (0) = 1 & , g^{e} = 1 + \frac{e}{1!} (g - 1) + \frac{e (e - 1)}{2!} (g - 1)^{2} + \dots \end{aligned}

在复合逆的计算中我们也会用到求幂函数．

Kronecker 代换

在分析时间复杂度之前我们先考虑如何作二元多项式乘法，一种想法是将系数「打包」，这一方法由 Kronecker 在 1882 年通过 $y \mapsto x^{N}$ 将 $R [x, y]$ 上的乘法缩减为 $R [x]$ 上的乘法，但是要求 $N$ 足够大．

不妨设 $\deg_{x} (A B) < N$ ，那么我们计算 $A (x, x^{N}) B (x, x^{N})$ 之后仍然可以还原出 $A (x, y) B (x, y)$ 且「打包」和「拆包」的时间为线性．

我们使用 Kronecker 代换再计算一元多项式乘法即可，不难发现在 $n$ 为二的幂时上述算法可以在 $O (M (n) \log n)$ 时间完成，因为每一次递归中 $y$ 的次数翻倍，但是 $x$ 的次数减半．

模板（P5373【模板】多项式复合函数）

代码相对于原算法作了一些简化及修改，使得代码更短．

#include <algorithm>
#include <cassert>
#include <cstdio>
#include <tuple>
#include <utility>
#include <vector>

using uint = unsigned;
using ull = unsigned long long;
constexpr uint MOD = 998244353;

constexpr uint PowMod(uint a, ull e) {
  for (uint res = 1;; a = (ull)a * a % MOD) {
    if (e & 1) res = (ull)res * a % MOD;
    if ((e /= 2) == 0) return res;
  }
}

constexpr uint InvMod(uint a) { return PowMod(a, MOD - 2); }

constexpr uint QUAD_NONRESIDUE = 3;
constexpr int LOG2_ORD = 23;  // __builtin_ctz(MOD - 1)
constexpr uint ZETA = PowMod(QUAD_NONRESIDUE, (MOD - 1) >> LOG2_ORD);
constexpr uint INV_ZETA = InvMod(ZETA);

std::vector<uint> GetFFTRoot(int n) {
  assert((n & (n - 1)) == 0);
  if (n / 2 == 0) return {};
  std::vector<uint> root(n / 2);
  root[0] = 1;
  for (int i = 0; (1 << i) < n / 2; ++i)
    root[1 << i] = PowMod(ZETA, 1LL << (LOG2_ORD - i - 2));
  for (int i = 1; i < n / 2; ++i)
    root[i] = (ull)root[i - (i & (i - 1))] * root[i & (i - 1)] % MOD;
  return root;
}

void Butterfly(uint a[], int n, const uint root[]) {
  assert((n & (n - 1)) == 0);
  for (int i = n; i >= 2; i /= 2)
    for (int j = 0; j < n; j += i)
      for (int k = j; k < j + i / 2; ++k) {
        const uint u = a[k];
        a[k + i / 2] = (ull)a[k + i / 2] * root[j / i] % MOD;
        if ((a[k] += a[k + i / 2]) >= MOD) a[k] -= MOD;
        if ((a[k + i / 2] = u + MOD - a[k + i / 2]) >= MOD) a[k + i / 2] -= MOD;
      }
}

void InvButterfly(uint a[], int n, const uint root[]) {
  assert((n & (n - 1)) == 0);
  for (int i = 2; i <= n; i *= 2)
    for (int j = 0; j < n; j += i)
      for (int k = j; k < j + i / 2; ++k) {
        const uint u = a[k];
        if ((a[k] += a[k + i / 2]) >= MOD) a[k] -= MOD;
        a[k + i / 2] = (ull)(u + MOD - a[k + i / 2]) * root[j / i] % MOD;
      }
}

void FFT(uint a[], int n, const uint root[]) { Butterfly(a, n, root); }

void InvFFT(uint a[], int n, const uint root[]) {
  InvButterfly(a, n, root);
  std::reverse(a + 1, a + n);
  const uint invn = InvMod(n);
  for (int i = 0; i < n; ++i) a[i] = (ull)a[i] * invn % MOD;
}

void TransposedFFT(uint a[], int n, const uint root[]) {
  InvButterfly(a, n, root);
}

void TransposedInvFFT(uint a[], int n, const uint root[]) {
  const uint invn = InvMod(n);
  for (int i = 0; i < n; ++i) a[i] = (ull)a[i] * invn % MOD;
  std::reverse(a + 1, a + n);
  Butterfly(a, n, root);
}

// 形式幂级数复合，求出 f(g) mod x^n 要求 g(0) = 0
std::vector<uint> FPSComp(std::vector<uint> f, std::vector<uint> g, int n) {
  assert(g.empty() || g[0] == 0);
  int N = 1;
  while (N < n) N *= 2;
  const std::vector<uint> root = GetFFTRoot(N * 4);
  // [y^(-1)] (f(y) / (-g(x) + y)) mod x^n in R[x]((y^(-1)))
  const auto KinoshitaLi = [&N, &root](auto&& KinoshitaLi, std::vector<uint>& P,
                                       std::vector<uint> Q, int d, int n) {
    assert((int)P.size() == N * 2);
    assert((int)Q.size() == N * 2);
    if (n == 1) return;
    Q.resize(N * 4);
    Q[N * 2] = 1;
    FFT(Q.data(), N * 4, root.data());
    if (n > 2) {
      std::vector<uint> V(N * 2);
      for (int i = 0; i < N * 4; i += 2) V[i / 2] = (ull)Q[i] * Q[i + 1] % MOD;
      InvFFT(V.data(), N * 2, root.data());
      assert(V[0] == 1);
      V[0] = 0;
      for (int i = 0; i < d * 2; ++i)
        fill_n(begin(V) + (i * n + n / 2), n / 2, 0u);
      KinoshitaLi(KinoshitaLi, P, std::move(V), d * 2, n / 2);
    }
    FFT(P.data(), N * 2, root.data());
    for (int i = 0; i < N * 4; i += 2)
      std::swap(Q[i] = (ull)Q[i] * P[i / 2] % MOD,
                Q[i + 1] = (ull)Q[i + 1] * P[i / 2] % MOD);
    InvFFT(Q.data(), N * 4, root.data());
    for (int i = 0; i < d; ++i)
      fill_n(copy_n(begin(Q) + (i + d) * n * 2, n, begin(P) + i * n * 2), n,
             0u);
  };
  f.resize(N * 2);
  g.resize(N * 2);
  for (int i = N - 1; i >= 0; --i) f[i * 2] = f[i], f[i * 2 + 1] = 0;
  for (int i = 0; i < N; ++i) g[i] = (g[i] != 0 ? MOD - g[i] : 0);
  fill_n(begin(g) + N, N, 0u);
  KinoshitaLi(KinoshitaLi, f, std::move(g), 1, N);
  f.resize(n);
  return f;
}

std::pair<std::vector<uint>, std::vector<uint>> GetFactorial(int n) {
  if (n == 0) return {};
  std::vector<uint> factorial(n), inv_factorial(n);
  factorial[0] = inv_factorial[0] = 1;
  if (n == 1) return {factorial, inv_factorial};
  std::vector<uint> inv(n);
  inv[1] = 1;
  for (int i = 2; i < n; ++i)
    inv[i] = (ull)(MOD - MOD / i) * inv[MOD % i] % MOD;
  for (int i = 1; i < n; ++i)
    factorial[i] = (ull)factorial[i - 1] * i % MOD,
    inv_factorial[i] = (ull)inv_factorial[i - 1] * inv[i] % MOD;
  return {factorial, inv_factorial};
}

// f(x) |-> f(x + c)
std::vector<uint> TaylorShift(std::vector<uint> f, uint c) {
  if (f.empty() || c == 0) return f;
  const int n = (int)f.size();
  std::vector<uint> factorial, inv_factorial;
  tie(factorial, inv_factorial) = GetFactorial(n);
  for (int i = 0; i < n; ++i) f[i] = (ull)f[i] * factorial[i] % MOD;
  std::vector<uint> g(n);
  uint cp = 1;
  for (int i = 0; i < n; ++i)
    g[i] = (ull)cp * inv_factorial[i] % MOD, cp = (ull)cp * c % MOD;
  int N = 1;
  while (N < n * 2 - 1) N *= 2;
  const std::vector<uint> root = GetFFTRoot(N);
  f.resize(N);
  g.resize(N);
  TransposedInvFFT(f.data(), N, root.data());
  FFT(g.data(), N, root.data());
  for (int i = 0; i < N; ++i) f[i] = (ull)f[i] * g[i] % MOD;
  TransposedFFT(f.data(), N, root.data());
  f.resize(n);
  for (int i = 0; i < n; ++i) f[i] = (ull)f[i] * inv_factorial[i] % MOD;
  return f;
}

// 多项式复合，求出 f(g) mod x^n
std::vector<uint> PolyComp(const std::vector<uint>& f,
                           const std::vector<uint>& g, int n) {
  if (g.empty() || g[0] == 0) return FPSComp(f, g, n);
  std::vector<uint> gg = g;
  gg[0] = 0;
  return FPSComp(TaylorShift(f, g[0]), std::move(gg), n);
}

int main() {
  int n, m;
  std::scanf("%d%d", &n, &m);
  std::vector<uint> f(n + 1), g(m + 1);
  for (int i = 0; i <= n; ++i) std::scanf("%u", &f[i]);
  for (int i = 0; i <= m; ++i) std::scanf("%u", &g[i]);
  const std::vector<uint> res = PolyComp(f, g, n + 1);
  for (int i = 0; i <= n; ++i) {
    if (i != 0) std::printf(" ");
    std::printf("%u", res[i]);
  }
  std::puts("");
  return 0;
}

形式幂级数的复合逆

现给出 $f \in x C [[x]]$ 且 $f^{'} (0) \neq 0$ ，求出 $g (x) mod x^{n}$ 满足 $f (g) \equiv g (f) \equiv x (\mod x^{n})$ ．

根据 Lagrange 反演，对于 $n > 1, k \geq 0$ 我们有

[x^{n - 1}] f (x)^{k} = \frac{k}{n - 1} [x^{n - 1 - k}] {(\frac{g (x)}{x})}^{- (n - 1)}

也就是我们如果能对 $k = 0, 1, \dots, n - 1$ 求出 $[x^{n - 1}] f (x)^{k}$ ，那么就可以求出其复合逆．

Kinoshita 和 Li 指出我们可以考虑二元有理函数

\frac{1}{1 - y \cdot f (x)} = \sum_{j \geq 0} f (x)^{j} y^{j}

且这个问题有一个更一般的形式即 Power Projection 问题：我们考虑计算

u := [x^{n - 1}] \frac{P (x, y)}{Q (x, y)} mod y^{m}

当 $n - 1 = 0$ 时显然有 $u = \frac{P (0, y)}{Q (0, y)} mod y^{m}$ ，否则我们有

\frac{P (x, y)}{Q (x, y)} = \frac{P (x, y) Q (- x, y)}{Q (x, y) Q (- x, y)} = \frac{U_{e} (x^{2}, y) + x U_{o} (x^{2}, y)}{V (x^{2}, y)}

那么

\begin{aligned} u & = {\begin{cases} [x^{n - 1}] \frac{U_{e} (x^{2}, y)}{V (x^{2}, y)} & if n - 1 is even, \\ [x^{n - 1}] \frac{x U_{o} (x^{2}, y)}{V (x^{2}, y)} & if n - 1 is odd. \end{cases} \\ = {\begin{cases} [x^{⌈ n / 2 ⌉ - 1}] \frac{U_{e} (x, y)}{V (x, y)} & if n - 1 is even, \\ [x^{⌈ n / 2 ⌉ - 1}] \frac{U_{o} (x, y)}{V (x, y)} & if n - 1 is odd. \end{cases} \end{aligned}

我们给出其伪代码：

\begin{array}{ll} Algorithm PowProj (P (x, y), Q (x, y), n, m) : \\ Input : P, Q \in C [x, y], n, m \in N_{> 0} . \\ Output : [x^{n - 1}] \frac{P (x, y)}{Q (x, y)} mod y^{m} . \\ Require : [x^{0} y^{0}] Q = 1 . \\ 1 & while n > 1 do \\ 2 & U (x, y) \leftarrow P (x, y) Q (- x, y) mod x^{n} mod y^{m} \\ 3 & if n - 1 is even then \\ 4 & P (x, y) \leftarrow \sum_{j = 0}^{⌈ n / 2 ⌉ + 1} ([x^{2 j}] U (x, y)) x^{j} \\ 5 & else \\ 6 & P (x, y) \leftarrow \sum_{j = 0}^{⌈ n / 2 ⌉ + 1} ([x^{2 j + 1}] U (x, y)) x^{j} \\ 7 & end if \\ 8 & V (x, y) \leftarrow Q (x, y) Q (- x, y) mod x^{n} mod y^{m} \\ 9 & Q (x, y) \leftarrow \sum_{j = 0}^{⌈ n / 2 ⌉ + 1} ([x^{2 j}] V (x, y)) x^{j} \\ 10 & n \leftarrow ⌈ n / 2 ⌉ \\ 11 & end while \\ 12 & return (\frac{P (0, y)}{Q (0, y)} mod y^{m}) \end{array}

同样的我们也可以直接导出 $Q (0, y)^{- 1}$ 而不需要计算形式幂级数的乘法逆元，那么复合逆的算法就是

\begin{array}{ll} Algorithm Rev (f (x), n) : \\ Input : f \in x C [[x]], f^{'} (0) \neq 0, n \in N_{\geq 2} . \\ Output : g (x) mod x^{n} such that f (g) \equiv g (f) \equiv x (\mod x^{n}) . \\ 1 & t \leftarrow f^{'} (0) \\ 2 & F (x) \leftarrow f (t^{- 1} x) \\ 3 & \sum_{k = 0}^{n - 1} a_{k} y^{k} \leftarrow PowProj (1, 1 - y \cdot F (x), n, n) \\ 4 & G (x) \leftarrow \sum_{k = 1}^{n - 1} \frac{n - 1}{k} a_{k} x^{n - 1 - k} \\ 5 & H (x) \leftarrow {(G (x)^{1 / (n - 1)})}^{- 1} mod x^{n - 1} \\ 6 & return ((t^{- 1} x) \circ (x \cdot H)) \end{array}

模板（P5809【模板】多项式复合逆）

代码相对于原算法作了一些简化及修改，使得代码更短．

#include <algorithm>
#include <cassert>
#include <cstdio>
#include <tuple>
#include <utility>
#include <vector>

using uint = unsigned;
using ull = unsigned long long;
constexpr uint MOD = 998244353;

constexpr uint PowMod(uint a, ull e) {
  for (uint res = 1;; a = (ull)a * a % MOD) {
    if (e & 1) res = (ull)res * a % MOD;
    if ((e /= 2) == 0) return res;
  }
}

constexpr uint InvMod(uint a) { return PowMod(a, MOD - 2); }

constexpr uint QUAD_NONRESIDUE = 3;
constexpr int LOG2_ORD = 23;  // __builtin_ctz(MOD - 1)
constexpr uint ZETA = PowMod(QUAD_NONRESIDUE, (MOD - 1) >> LOG2_ORD);
constexpr uint INV_ZETA = InvMod(ZETA);

std::vector<uint> GetFFTRoot(int n) {
  assert((n & (n - 1)) == 0);
  if (n / 2 == 0) return {};
  std::vector<uint> root(n / 2);
  root[0] = 1;
  for (int i = 0; (1 << i) < n / 2; ++i)
    root[1 << i] = PowMod(ZETA, 1LL << (LOG2_ORD - i - 2));
  for (int i = 1; i < n / 2; ++i)
    root[i] = (ull)root[i - (i & (i - 1))] * root[i & (i - 1)] % MOD;
  return root;
}

void Butterfly(uint a[], int n, const uint root[]) {
  assert((n & (n - 1)) == 0);
  for (int i = n; i >= 2; i /= 2)
    for (int j = 0; j < n; j += i)
      for (int k = j; k < j + i / 2; ++k) {
        const uint u = a[k];
        a[k + i / 2] = (ull)a[k + i / 2] * root[j / i] % MOD;
        if ((a[k] += a[k + i / 2]) >= MOD) a[k] -= MOD;
        if ((a[k + i / 2] = u + MOD - a[k + i / 2]) >= MOD) a[k + i / 2] -= MOD;
      }
}

void InvButterfly(uint a[], int n, const uint root[]) {
  assert((n & (n - 1)) == 0);
  for (int i = 2; i <= n; i *= 2)
    for (int j = 0; j < n; j += i)
      for (int k = j; k < j + i / 2; ++k) {
        const uint u = a[k];
        if ((a[k] += a[k + i / 2]) >= MOD) a[k] -= MOD;
        a[k + i / 2] = (ull)(u + MOD - a[k + i / 2]) * root[j / i] % MOD;
      }
}

void FFT(uint a[], int n, const uint root[]) { Butterfly(a, n, root); }

void InvFFT(uint a[], int n, const uint root[]) {
  InvButterfly(a, n, root);
  std::reverse(a + 1, a + n);
  const uint invn = InvMod(n);
  for (int i = 0; i < n; ++i) a[i] = (ull)a[i] * invn % MOD;
}

void TransposedFFT(uint a[], int n, const uint root[]) {
  InvButterfly(a, n, root);
}

void TransposedInvFFT(uint a[], int n, const uint root[]) {
  const uint invn = InvMod(n);
  for (int i = 0; i < n; ++i) a[i] = (ull)a[i] * invn % MOD;
  std::reverse(a + 1, a + n);
  Butterfly(a, n, root);
}

// 形式幂级数复合，求出 f(g) mod x^n 要求 g(0) = 0
std::vector<uint> FPSComp(std::vector<uint> f, std::vector<uint> g, int n) {
  assert(g.empty() || g[0] == 0);
  int N = 1;
  while (N < n) N *= 2;
  const std::vector<uint> root = GetFFTRoot(N * 4);
  // [y^(-1)] (f(y) / (-g(x) + y)) mod x^n in R[x]((y^(-1)))
  const auto KinoshitaLi = [&N, &root](auto&& KinoshitaLi, std::vector<uint>& P,
                                       std::vector<uint> Q, int d, int n) {
    assert((int)P.size() == N * 2);
    assert((int)Q.size() == N * 2);
    if (n == 1) return;
    Q.resize(N * 4);
    Q[N * 2] = 1;
    FFT(Q.data(), N * 4, root.data());
    if (n > 2) {
      std::vector<uint> V(N * 2);
      for (int i = 0; i < N * 4; i += 2) V[i / 2] = (ull)Q[i] * Q[i + 1] % MOD;
      InvFFT(V.data(), N * 2, root.data());
      assert(V[0] == 1);
      V[0] = 0;
      for (int i = 0; i < d * 2; ++i)
        fill_n(begin(V) + (i * n + n / 2), n / 2, 0u);
      KinoshitaLi(KinoshitaLi, P, std::move(V), d * 2, n / 2);
    }
    FFT(P.data(), N * 2, root.data());
    for (int i = 0; i < N * 4; i += 2)
      std::swap(Q[i] = (ull)Q[i] * P[i / 2] % MOD,
                Q[i + 1] = (ull)Q[i + 1] * P[i / 2] % MOD);
    InvFFT(Q.data(), N * 4, root.data());
    for (int i = 0; i < d; ++i)
      fill_n(copy_n(begin(Q) + (i + d) * n * 2, n, begin(P) + i * n * 2), n,
             0u);
  };
  f.resize(N * 2);
  g.resize(N * 2);
  for (int i = N - 1; i >= 0; --i) f[i * 2] = f[i], f[i * 2 + 1] = 0;
  for (int i = 0; i < N; ++i) g[i] = (g[i] != 0 ? MOD - g[i] : 0);
  fill_n(begin(g) + N, N, 0u);
  KinoshitaLi(KinoshitaLi, f, std::move(g), 1, N);
  f.resize(n);
  return f;
}

std::pair<std::vector<uint>, std::vector<uint>> GetFactorial(int n) {
  if (n == 0) return {};
  std::vector<uint> factorial(n), inv_factorial(n);
  factorial[0] = inv_factorial[0] = 1;
  if (n == 1) return {factorial, inv_factorial};
  std::vector<uint> inv(n);
  inv[1] = 1;
  for (int i = 2; i < n; ++i)
    inv[i] = (ull)(MOD - MOD / i) * inv[MOD % i] % MOD;
  for (int i = 1; i < n; ++i)
    factorial[i] = (ull)factorial[i - 1] * i % MOD,
    inv_factorial[i] = (ull)inv_factorial[i - 1] * inv[i] % MOD;
  return {factorial, inv_factorial};
}

// 二项卷积
// Kinoshita--Li 算法作者之一的李白天在 2021 年的
// 集训队论文《信息学竞赛中的生成函数计算理论框架》中提到
// Taylor Shift 算法为二项卷积的转置
std::vector<uint> BinConv(std::vector<uint> f, uint c) {
  if (f.empty() || c == 0) return f;
  const int n = (int)f.size();
  std::vector<uint> factorial, inv_factorial;
  tie(factorial, inv_factorial) = GetFactorial(n);
  for (int i = 0; i < n; ++i) f[i] = (ull)f[i] * inv_factorial[i] % MOD;
  std::vector<uint> g(n);
  uint cp = 1;
  for (int i = 0; i < n; ++i)
    g[i] = (ull)cp * inv_factorial[i] % MOD, cp = (ull)cp * c % MOD;
  int N = 1;
  while (N < n * 2 - 1) N *= 2;
  const std::vector<uint> root = GetFFTRoot(N);
  f.resize(N);
  g.resize(N);
  FFT(f.data(), N, root.data());
  FFT(g.data(), N, root.data());
  for (int i = 0; i < N; ++i) f[i] = (ull)f[i] * g[i] % MOD;
  InvFFT(f.data(), N, root.data());
  f.resize(n);
  for (int i = 0; i < n; ++i) f[i] = (ull)f[i] * factorial[i] % MOD;
  return f;
}

// Power Projection: [x^(n-1)] (f g^i) for i = 0, ..., n-1
// 不难发现该算法就是 PolyComp 算法的转置
std::vector<uint> PowProj(std::vector<uint> f, std::vector<uint> g, int n) {
  const uint c = g.empty() ? 0u : g[0];
  if (!g.empty()) g[0] = 0;
  int N = 1;
  while (N < n) N *= 2;
  const std::vector<uint> root = GetFFTRoot(N * 4);
  // [x^(n-1)] (f(x) / (-g(x) + y)) in R[x]((y^(-1)))
  const auto KinoshitaLi = [&N, &root](std::vector<uint>& P,
                                       std::vector<uint>& Q, int d, int n) {
    assert((int)P.size() == N * 2);
    assert((int)Q.size() == N * 2);
    P.insert(begin(P), N * 2, 0u);
    uint* const nP = P.data() + N * 2;
    Q.resize(N * 4);
    for (; n > 1; d *= 2, n /= 2) {
      Q[N * 2] = 1;
      TransposedInvFFT(P.data(), N * 4, root.data());
      FFT(Q.data(), N * 4, root.data());
      for (int i = N * 4 - 2; i >= 0; i -= 2)
        nP[i / 2] = ((ull)P[i] * Q[i + 1] + (ull)P[i + 1] * Q[i]) % MOD;
      TransposedFFT(nP, N * 2, root.data());
      if (n == 2) break;
      for (int i = 0; i < N * 4; i += 2) Q[i / 2] = (ull)Q[i] * Q[i + 1] % MOD;
      InvFFT(Q.data(), N * 2, root.data());
      assert(Q[0] == 1);
      Q[0] = 0;
      for (int i = 0; i < d * 2; ++i) {
        std::fill_n(nP + i * n, n / 2, 0u);
        fill_n(begin(Q) + (i * n + n / 2), n / 2, 0u);
      }
      fill_n(begin(P), N * 2, 0u);
      fill_n(begin(Q) + N * 2, N * 2, 0u);
    }
    P.erase(begin(P), begin(P) + N * 2);
  };
  f.insert(begin(f), N - n, 0u);
  f.resize(N);
  reverse(begin(f), end(f));
  f.insert(begin(f), N, 0u);
  g.resize(N * 2);
  for (int i = 0; i < N; ++i) g[i] = (g[i] != 0 ? MOD - g[i] : 0);
  fill_n(begin(g) + N, N, 0u);
  KinoshitaLi(f, g, 1, N);
  for (int i = 0; i < N; ++i) f[i] = f[i * 2 + 1];
  f.resize(n);
  return BinConv(std::move(f), c);
}

// 形式幂级数幂函数，计算 g^e mod x^n 要求 g(0) = 1
std::vector<uint> FPSPow1(std::vector<uint> g, uint e, int n) {
  assert(!g.empty() && g[0] == 1);
  if (n == 1) return std::vector<uint>{1u};
  std::vector<uint> inv(n), f(n);
  inv[1] = f[0] = 1;
  for (int i = 2; i < n; ++i)
    inv[i] = (ull)(MOD - MOD / i) * inv[MOD % i] % MOD;
  for (int i = 1; i < n; ++i)
    f[i] = (ull)f[i - 1] * (e + MOD + 1 - i) % MOD * inv[i] % MOD;
  g[0] = 0;
  return FPSComp(std::move(f), std::move(g), n);
}

std::vector<uint> FPSRev(std::vector<uint> f, int n) {
  assert(f.size() >= 2 && f[0] == 0 && f[1] != 0);
  if (n == 1) return std::vector<uint>{0u};
  f.resize(n);
  const uint invf1 = InvMod(f[1]);
  uint invf1p = 1;
  for (int i = 0; i < n; ++i)
    f[i] = (ull)f[i] * invf1p % MOD, invf1p = (ull)invf1p * invf1 % MOD;
  std::vector<uint> inv(n);
  inv[1] = 1;
  for (int i = 2; i < n; ++i)
    inv[i] = (ull)(MOD - MOD / i) * inv[MOD % i] % MOD;
  std::vector<uint> proj = PowProj(std::vector<uint>{1u}, std::move(f), n);
  for (int i = 1; i < n; ++i)
    proj[i] = (ull)proj[i] * (n - 1) % MOD * inv[i] % MOD;
  reverse(begin(proj), end(proj));
  std::vector<uint> res = FPSPow1(proj, MOD - inv[n - 1], n - 1);
  for (int i = 0; i < n - 1; ++i) res[i] = (ull)res[i] * invf1 % MOD;
  res.insert(begin(res), 0);
  return res;
}

int main() {
  int n;
  std::scanf("%d", &n);
  std::vector<uint> f(n);
  for (int i = 0; i < n; ++i) std::scanf("%u", &f[i]);
  const std::vector<uint> res = FPSRev(f, n);
  for (int i = 0; i < n; ++i) {
    if (i != 0) std::printf(" ");
    std::printf("%u", res[i]);
  }
  std::puts("");
  return 0;
}

由转置原理导出

Power Projection 问题是 Modular Composition 的转置，Kinoshita 和 Li 指出我们前文的复合算法可以由 Power Projection 算法直接转置得到．同样的，如果优化可以应用于 Power Projection 算法，其也可以应用于 Modular Composition 算法．我们省略细节．

参考文献

Yasunori Kinoshita, Baitian Li.Power Series Composition in Near-Linear Time. FOCS 2024.
Alin Bostan, Ryuhei Mori.A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence. SOSA 2021: 118-132
R. P. Brent and H. T. Kung. 1978.Fast Algorithms for Manipulating Formal Power Series. J. ACM 25, 4 (Oct. 1978), 581–595.
Daniel J. Bernstein. "Fast multiplication and its applications." Pages 325–384 in Algorithmic number theory: lattices, number fields, curves and cryptography, edited by Joe Buhler, Peter Stevenhagen, Cambridge University Press, 2008, ISBN 978-0521808545.

本页面最近更新：2026/6/21 16:39:02，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：hly1204, Tiphereth-A, c-forrest
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用