跳转至

# 形式幂级数复合|复合逆

形式幂级数的复合和复合逆也是常见的形式幂级数操作，对于没有特殊性质的 𝑓 f 之前我们一直使用的多是 𝑂(𝑛 2)O(n 2) 的算法来计算 𝑓(𝑔)m o d 𝑥 𝑛 f(g)mod x n 其中 𝑓 ∈ℂ[[𝑥]],𝑔 ∈𝑥 ℂ[[𝑥]]f∈C[[x]],g∈x C[[x]]，但是因为效率较低应用较少．我们介绍 Kinoshita–Li 的 𝑂(𝖬(𝑛)l o g⁡𝑛)O(M(n)log⁡n) 的算法，其中 𝑂(𝖬(𝑛))O(M(n)) 为两个次数为 𝑂(𝑛)O(n) 的多项式相乘的时间．

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

若要计算 𝑓(𝑔(𝑥))m o d 𝑥 𝑛 f(g(x))mod x n 那么需要 𝑓(𝑔(𝑥))f(g(x)) 的每一项系数都是有限项之和，所以之前要求 𝑓(𝑥) ∈ℂ[[𝑥]],𝑔(𝑥) ∈𝑥 ℂ[[𝑥]]f(x)∈C[[x]],g(x)∈x C[[x]]，而如果 𝑓(𝑥),𝑔(𝑥) ∈ℂ[𝑥]f(x),g(x)∈C[x] 也可以满足这个条件．因为我们需要将 𝑓(𝑔(𝑥))f(g(x)) 的系数截断，不妨直接考虑 𝑓(𝑥),𝑔(𝑥)f(x),g(x) 都是多项式的情况．对于 𝑓(𝑥) =∑𝑛−1 𝑗=0 𝑓 𝑗 𝑥 𝑗 f(x)=∑j=0 n−1 f j x j，有

𝑓(𝑔(𝑥))=𝑛−1∑𝑗=0 𝑓 𝑗 𝑔(𝑥)𝑗 f(g(x))=∑j=0 n−1 f j g(x)j

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

𝑓(𝑦−1)1−𝑦⋅𝑔(𝑥)=∑𝑗≥0(⋯+𝑓 𝑗 𝑦−𝑗+⋯)𝑔(𝑥)𝑗 𝑦 𝑗 𝑓(𝑔(𝑥))=[𝑦 0]𝑓(𝑦−1)1−𝑦⋅𝑔(𝑥)f(y−1)1−y⋅g(x)=∑j≥0(⋯+f j y−j+⋯)g(x)j y j f(g(x))=[y 0]f(y−1)1−y⋅g(x)

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

𝑃(𝑦)𝑄(𝑥,𝑦)m o d 𝑥 𝑛=(𝑃(𝑦)𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛=(𝑃(𝑦)𝑉(𝑥 2,𝑦)m o d 𝑥 𝑛)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛=(𝑃(𝑦)𝑉(𝑧,𝑦)m o d 𝑧⌈𝑛/2⌉)∣𝑧=𝑥 2 𝑄(−𝑥,𝑦)m o d 𝑥 𝑛 P(y)Q(x,y)mod x n=(P(y)Q(x,y)Q(−x,y)mod x n)Q(−x,y)mod x n=(P(y)V(x 2,y)mod x n)Q(−x,y)mod x n=(P(y)V(z,y)mod z⌈n/2⌉)|z=x 2 Q(−x,y)mod x n

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

𝑃(𝑦)𝑄(𝑥,𝑦)m o d 𝑥=𝑃(𝑦)𝑄(0,𝑦)∈ℂ((𝑦))P(y)Q(x,y)mod x=P(y)Q(0,y)∈C((y))

在计算 𝑃(𝑦)𝑉(𝑧,𝑦)m o d 𝑧⌈𝑛/2⌉ ∈ℂ[𝑧]((𝑦))P(y)V(z,y)mod z⌈n/2⌉∈C[z]((y)) 时我们不需要保留所有 𝑦 y 的系数，因为最后我们只需要提取 𝑦 0 y 0 的系数，所以 𝑦>0 y>0 的系数是不需要的，而因为求出前者之后需要将其乘以若干个形如 𝑄( −𝑥,𝑦) ∈ℂ[𝑥,𝑦]Q(−x,y)∈C[x,y] 的「多项式」，所以只需要保留对于 𝑦 0 y 0 有贡献的系数即可．我们准备好给出伪代码：

𝐀 𝐥 𝐠 𝐨 𝐫 𝐢 𝐭 𝐡 𝐦 𝖢 𝗈 𝗆 𝗉⁡(𝑃(𝑦),𝑄(𝑥,𝑦),𝑛,𝑚):𝐈 𝐧 𝐩 𝐮 𝐭: 𝑃=∑0≤𝑗<𝑛 𝑝 𝑗 𝑦−𝑗∈ℂ((𝑦)),𝑄∈ℂ[𝑥,𝑦],𝑛,𝑚∈ℕ>0.𝐎 𝐮 𝐭 𝐩 𝐮 𝐭: [𝑦(−𝑚,0]]𝑃(𝑦)𝑄(𝑥,𝑦)m o d 𝑥 𝑛.𝐑 𝐞 𝐪 𝐮 𝐢 𝐫 𝐞: [𝑥 0 𝑦 0]𝑄=1.1 𝐢 𝐟 𝑛=1 𝐭 𝐡 𝐞 𝐧 𝐫 𝐞 𝐭 𝐮 𝐫 𝐧 ([𝑦−𝑚+1]𝑃(𝑦)𝑄(0,𝑦),…,[𝑦 0]𝑃(𝑦)𝑄(0,𝑦))2 𝑉(𝑥 2,𝑦)←𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛 3 𝑑←d e g 𝑦⁡𝑄(−𝑥,𝑦)4(𝑡−(𝑚+𝑑)+1,…,𝑡 0)←𝖢 𝗈 𝗆 𝗉⁡(𝑃(𝑦),𝑉(𝑥,𝑦),⌈𝑛/2⌉,𝑚+𝑑)5 𝑇(𝑥,𝑦)←∑0 𝑗=−(𝑚+𝑑)+1 𝑡 𝑗 𝑦 𝑗 6 𝑈(𝑥,𝑦)=∑𝑑 𝑗=−(𝑚+𝑑)+1 𝑢 𝑗 𝑦 𝑗←𝑇(𝑥 2,𝑦)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛 7 𝐫 𝐞 𝐭 𝐮 𝐫 𝐧 (𝑢−𝑚+1,…,𝑢 0)Algorithm Comp⁡(P(y),Q(x,y),n,m):Input: P=∑0≤j<n p j y−j∈C((y)),Q∈C[x,y],n,m∈N>0.Output: [y(−m,0]]P(y)Q(x,y)mod x n.Require: [x 0 y 0]Q=1.1 if n=1 then return ([y−m+1]P(y)Q(0,y),…,[y 0]P(y)Q(0,y))2 V(x 2,y)←Q(x,y)Q(−x,y)mod x n 3 d←deg y⁡Q(−x,y)4(t−(m+d)+1,…,t 0)←Comp⁡(P(y),V(x,y),⌈n/2⌉,m+d)5 T(x,y)←∑j=−(m+d)+1 0 t j y j 6 U(x,y)=∑j=−(m+d)+1 d u j y j←T(x 2,y)Q(−x,y)mod x n 7 return (u−m+1,…,u 0)

那么我们有

𝑓(𝑔(𝑥))m o d 𝑥 𝑛=𝖢 𝗈 𝗆 𝗉⁡(𝑓(𝑦−1),1−𝑦⋅𝑔(𝑥),m a x{1+d e g⁡𝑓,𝑛},1)m o d 𝑥 𝑛 f(g(x))mod x n=Comp⁡(f(y−1),1−y⋅g(x),max{1+deg⁡f,n},1)mod x n

注意第三个参数是因为 𝑔(0)g(0) 可能不为零，如果 d e g⁡𝑓 ≥𝑛 deg⁡f≥n 此时不能截断 𝑓(𝑥)f(x) 来计算 𝑓(𝑔(𝑥))f(g(x))，我们也可以选择计算 𝑓(𝑔) =𝑓 ∘(𝑥+𝑔(0)) ∘(𝑔−𝑔(0))f(g)=f∘(x+g(0))∘(g−g(0))，此时可以取 𝐹 :=𝑓(𝑥+𝑔(0))m o d 𝑥 𝑛 F:=f(x+g(0))mod x n 和 𝐺 :=𝑔 −𝑔(0)G:=g−g(0) 转而计算 𝖢 𝗈 𝗆 𝗉⁡(𝐹(𝑦−1),1−𝑦⋅𝐺(𝑥),𝑛,1)Comp⁡(F(y−1),1−y⋅G(x),n,1)．

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

## 常见的特殊形式复合

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

𝑔(0)=1, 𝑔−1=1+(1−𝑔)+(1−𝑔)2+⋯𝑔(0)=1, l o g⁡𝑔=−1−𝑔 1−(1−𝑔)2 2−(1−𝑔)3 3−⋯𝑔(0)=0, e x p⁡𝑔=1+𝑔 1!+𝑔 2 2!+𝑔 3 3!+⋯𝑔(0)=1, 𝑔 𝑒=1+𝑒 1(𝑔−1)+𝑒(𝑒−1)2(𝑔−1)2+⋯g(0)=1, g−1=1+(1−g)+(1−g)2+⋯g(0)=1, log⁡g=−1−g 1−(1−g)2 2−(1−g)3 3−⋯g(0)=0, exp⁡g=1+g 1!+g 2 2!+g 3 3!+⋯g(0)=1, g e=1+e 1(g−1)+e(e−1)2(g−1)2+⋯

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

### Kronecker 代换

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

不妨设 d e g 𝑥⁡(𝐴 𝐵) <𝑁 deg x⁡(A B)<N，那么我们计算 𝐴(𝑥,𝑥 𝑁)𝐵(𝑥,𝑥 𝑁)A(x,x N)B(x,x N) 之后仍然可以还原出 𝐴(𝑥,𝑦)𝐵(𝑥,𝑦)A(x,y)B(x,y) 且「打包」和「拆包」的时间为线性．

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

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

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

      1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
     89
     90
     91
     92
     93
     94
     95
     96
     97
     98
     99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    179
    180
    181
    182
    183
    #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);

    // 返回做 n 长 FFT 所需的单位根数组，长度为一半
    std::pair<std::vector<uint>, std::vector<uint>> GetFFTRoot(int n) {
      assert((n & (n - 1)) == 0);
      if (n / 2 == 0) return {};
      std::vector<uint> root(n / 2), inv_root(n / 2);
      root[0] = inv_root[0] = 1;
      for (int i = 0; (1 << i) < n / 2; ++i)
        root[1 << i] = PowMod(ZETA, 1LL << (LOG2_ORD - i - 2)),
                  inv_root[1 << i] = PowMod(INV_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,
        inv_root[i] =
            (ull)inv_root[i - (i & (i - 1))] * inv_root[i & (i - 1)] % MOD;
      return {root, inv_root};
    }

    void Butterfly(int n, uint a[], 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(int n, uint a[], 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(int n, uint a[], const uint root[]) { Butterfly(n, a, root); }

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

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

    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 = f.size();
      std::vector<uint> factorial, inv_factorial;
      std::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 len = 1;
      while (len < n * 2 - 1) len *= 2;
      std::vector<uint> root, inv_root;
      std::tie(root, inv_root) = GetFFTRoot(len);
      f.resize(len);
      g.resize(len);
      FFT(len, f.data(), inv_root.data());
      FFT(len, g.data(), root.data());
      for (int i = 0; i < len; ++i) f[i] = (ull)f[i] * g[i] % MOD;
      InvFFT(len, f.data(), 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> PolyComposition(const std::vector<uint> &f,
                                      const std::vector<uint> &g, int n) {
      if (g.empty() || g[0] == 0) return FPSComposition(f, g, n);
      std::vector<uint> gg = g;
      gg[0] = 0;
      return FPSComposition(TaylorShift(f, g[0]), 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 = PolyComposition(f, g, n + 1);
      for (int i = 0; i <= n; ++i) std::printf("%u%c", res[i], " \n"[i == n]);
      return 0;
    }

## 形式幂级数的复合逆

现给出 𝑓 ∈𝑥 ℂ[[𝑥]]f∈x C[[x]] 且 𝑓′(0) ≠0 f′(0)≠0，求出 𝑔(𝑥)m o d 𝑥 𝑛 g(x)mod x n 满足 𝑓(𝑔) ≡𝑔(𝑓) ≡𝑥(m o d 𝑥 𝑛)f(g)≡g(f)≡x(mod x n)．

根据 Lagrange 反演，对于 𝑛 >1,𝑘 ≥0 n>1,k≥0 我们有

[𝑥 𝑛−1]𝑓(𝑥)𝑘=𝑘 𝑛−1[𝑥 𝑛−1−𝑘](𝑔(𝑥)𝑥)−(𝑛−1)[x n−1]f(x)k=k n−1[x n−1−k](g(x)x)−(n−1)

也就是我们如果能对 𝑘 =0,1,…,𝑛 −1 k=0,1,…,n−1 求出 [𝑥 𝑛−1]𝑓(𝑥)𝑘[x n−1]f(x)k，那么就可以求出其复合逆．

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

1 1−𝑦⋅𝑓(𝑥)=∑𝑗≥0 𝑓(𝑥)𝑗 𝑦 𝑗 1 1−y⋅f(x)=∑j≥0 f(x)j y j

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

𝑢:=[𝑥 𝑛−1]𝑃(𝑥,𝑦)𝑄(𝑥,𝑦)m o d 𝑦 𝑚 u:=[x n−1]P(x,y)Q(x,y)mod y m

当 𝑛 −1 =0 n−1=0 时显然有 𝑢 =𝑃(0,𝑦)𝑄(0,𝑦)m o d 𝑦 𝑚 u=P(0,y)Q(0,y)mod y m，否则我们有

𝑃(𝑥,𝑦)𝑄(𝑥,𝑦)=𝑃(𝑥,𝑦)𝑄(−𝑥,𝑦)𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)=𝑈 𝑒(𝑥 2,𝑦)+𝑥 𝑈 𝑜(𝑥 2,𝑦)𝑉(𝑥 2,𝑦)P(x,y)Q(x,y)=P(x,y)Q(−x,y)Q(x,y)Q(−x,y)=U e(x 2,y)+x U o(x 2,y)V(x 2,y)

那么

𝑢=⎧{ { {⎨{ { {⎩[𝑥 𝑛−1]𝑈 𝑒(𝑥 2,𝑦)𝑉(𝑥 2,𝑦) i f 𝑛−1 i s e v e n,[𝑥 𝑛−1]𝑥 𝑈 𝑜(𝑥 2,𝑦)𝑉(𝑥 2,𝑦) i f 𝑛−1 i s o d d.=⎧{ { {⎨{ { {⎩[𝑥⌈𝑛/2⌉−1]𝑈 𝑒(𝑥,𝑦)𝑉(𝑥,𝑦) i f 𝑛−1 i s e v e n,[𝑥⌈𝑛/2⌉−1]𝑈 𝑜(𝑥,𝑦)𝑉(𝑥,𝑦) i f 𝑛−1 i s o d d.u={[x n−1]U e(x 2,y)V(x 2,y) if n−1 is even,[x n−1]x U o(x 2,y)V(x 2,y) if n−1 is odd.={[x⌈n/2⌉−1]U e(x,y)V(x,y) if n−1 is even,[x⌈n/2⌉−1]U o(x,y)V(x,y) if n−1 is odd.

我们给出其伪代码：

𝐀 𝐥 𝐠 𝐨 𝐫 𝐢 𝐭 𝐡 𝐦 𝖯 𝗈 𝗐 𝖯 𝗋 𝗈 𝗃⁡(𝑃(𝑥,𝑦),𝑄(𝑥,𝑦),𝑛,𝑚):𝐈 𝐧 𝐩 𝐮 𝐭: 𝑃,𝑄∈ℂ[𝑥,𝑦],𝑛,𝑚∈ℕ>0.𝐎 𝐮 𝐭 𝐩 𝐮 𝐭: [𝑥 𝑛−1]𝑃(𝑥,𝑦)𝑄(𝑥,𝑦)m o d 𝑦 𝑚.𝐑 𝐞 𝐪 𝐮 𝐢 𝐫 𝐞: [𝑥 0 𝑦 0]𝑄=1.1 𝐰 𝐡 𝐢 𝐥 𝐞 𝑛>1 𝐝 𝐨 2 𝑈(𝑥,𝑦)←𝑃(𝑥,𝑦)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛 m o d 𝑦 𝑚 3 𝐢 𝐟 𝑛−1 i s e v e n 𝐭 𝐡 𝐞 𝐧 4 𝑃(𝑥,𝑦)←∑⌈𝑛/2⌉+1 𝑗=0([𝑥 2 𝑗]𝑈(𝑥,𝑦))𝑥 𝑗 5 𝐞 𝐥 𝐬 𝐞 6 𝑃(𝑥,𝑦)←∑⌈𝑛/2⌉+1 𝑗=0([𝑥 2 𝑗+1]𝑈(𝑥,𝑦))𝑥 𝑗 7 𝐞 𝐧 𝐝 𝐢 𝐟 8 𝑉(𝑥,𝑦)←𝑄(𝑥,𝑦)𝑄(−𝑥,𝑦)m o d 𝑥 𝑛 m o d 𝑦 𝑚 9 𝑄(𝑥,𝑦)←∑⌈𝑛/2⌉+1 𝑗=0([𝑥 2 𝑗]𝑉(𝑥,𝑦))𝑥 𝑗 1 0 𝑛←⌈𝑛/2⌉1 1 𝐞 𝐧 𝐝 𝐰 𝐡 𝐢 𝐥 𝐞 1 2 𝐫 𝐞 𝐭 𝐮 𝐫 𝐧 (𝑃(0,𝑦)𝑄(0,𝑦)m o d 𝑦 𝑚)Algorithm PowProj⁡(P(x,y),Q(x,y),n,m):Input: P,Q∈C[x,y],n,m∈N>0.Output: [x n−1]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)←P(x,y)Q(−x,y)mod x n mod y m 3 if n−1 is even then 4 P(x,y)←∑j=0⌈n/2⌉+1([x 2 j]U(x,y))x j 5 else 6 P(x,y)←∑j=0⌈n/2⌉+1([x 2 j+1]U(x,y))x j 7 end if 8 V(x,y)←Q(x,y)Q(−x,y)mod x n mod y m 9 Q(x,y)←∑j=0⌈n/2⌉+1([x 2 j]V(x,y))x j 10 n←⌈n/2⌉11 end while 12 return (P(0,y)Q(0,y)mod y m)

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

𝐀 𝐥 𝐠 𝐨 𝐫 𝐢 𝐭 𝐡 𝐦 𝖱 𝖾 𝗏⁡(𝑓(𝑥),𝑛):𝐈 𝐧 𝐩 𝐮 𝐭: 𝑓∈𝑥 ℂ[[𝑥]],𝑓′(0)≠0,𝑛∈ℕ≥2.𝐎 𝐮 𝐭 𝐩 𝐮 𝐭: 𝑔(𝑥)m o d 𝑥 𝑛 s u c h t h a t 𝑓(𝑔)≡𝑔(𝑓)≡𝑥(m o d 𝑥 𝑛).1 𝑡←𝑓′(0)2 𝐹(𝑥)←𝑓(𝑡−1 𝑥)3∑𝑛−1 𝑘=0 𝑎 𝑘 𝑦 𝑘←𝖯 𝗈 𝗐 𝖯 𝗋 𝗈 𝗃⁡(1,1−𝑦⋅𝐹(𝑥),𝑛,𝑛)4 𝐺(𝑥)←∑𝑛−1 𝑘=1 𝑛−1 𝑘 𝑎 𝑘 𝑥 𝑛−1−𝑘 5 𝐻(𝑥)←(𝐺(𝑥)1/(𝑛−1))−1 m o d 𝑥 𝑛−1 6 𝐫 𝐞 𝐭 𝐮 𝐫 𝐧 ((𝑡−1 𝑥)∘(𝑥⋅𝐻))Algorithm Rev⁡(f(x),n):Input: f∈x C[[x]],f′(0)≠0,n∈N≥2.Output: g(x)mod x n such that f(g)≡g(f)≡x(mod x n).1 t←f′(0)2 F(x)←f(t−1 x)3∑k=0 n−1 a k y k←PowProj⁡(1,1−y⋅F(x),n,n)4 G(x)←∑k=1 n−1 n−1 k a k x n−1−k 5 H(x)←(G(x)1/(n−1))−1 mod x n−1 6 return ((t−1 x)∘(x⋅H))
模板（P5809【模板】多项式复合逆）

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

      1
      2
      3
      4
      5
      6
      7
      8
      9
     10
     11
     12
     13
     14
     15
     16
     17
     18
     19
     20
     21
     22
     23
     24
     25
     26
     27
     28
     29
     30
     31
     32
     33
     34
     35
     36
     37
     38
     39
     40
     41
     42
     43
     44
     45
     46
     47
     48
     49
     50
     51
     52
     53
     54
     55
     56
     57
     58
     59
     60
     61
     62
     63
     64
     65
     66
     67
     68
     69
     70
     71
     72
     73
     74
     75
     76
     77
     78
     79
     80
     81
     82
     83
     84
     85
     86
     87
     88
     89
     90
     91
     92
     93
     94
     95
     96
     97
     98
     99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    179
    180
    181
    182
    183
    184
    185
    186
    187
    188
    189
    190
    191
    192
    193
    194
    195
    196
    197
    198
    199
    200
    201
    202
    203
    204
    205
    206
    207
    208
    209
    210
    211
    212
    213
    214
    215
    216
    217
    218
    219
    220
    221
    222
    223
    224
    225
    #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);

    // 返回做 n 长 FFT 所需的单位根数组，长度为一半
    std::pair<std::vector<uint>, std::vector<uint>> GetFFTRoot(int n) {
      assert((n & (n - 1)) == 0);
      if (n / 2 == 0) return {};
      std::vector<uint> root(n / 2), inv_root(n / 2);
      root[0] = inv_root[0] = 1;
      for (int i = 0; (1 << i) < n / 2; ++i)
        root[1 << i] = PowMod(ZETA, 1LL << (LOG2_ORD - i - 2)),
                  inv_root[1 << i] = PowMod(INV_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,
        inv_root[i] =
            (ull)inv_root[i - (i & (i - 1))] * inv_root[i & (i - 1)] % MOD;
      return {root, inv_root};
    }

    void Butterfly(int n, uint a[], 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(int n, uint a[], 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(int n, uint a[], const uint root[]) { Butterfly(n, a, root); }

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

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

    // Power Projection: [x^(n-1)] (fg^i) for i=0,..,n-1 要求 g(0) = 0
    std::vector<uint> PowerProjection(std::vector<uint> f, std::vector<uint> g,
                                      int n) {
      assert(g.empty() || g[0] == 0);
      int len = 1;
      while (len < n) len *= 2;
      std::vector<uint> root, inv_root;
      std::tie(root, inv_root) = GetFFTRoot(len * 4);
      // [x^(n-1)] (f(x) / (-g(x) + y)) in R[x]((y^(-1)))
      const auto KinoshitaLi = [&](auto &&KinoshitaLi, std::vector<uint> P,
                                   std::vector<uint> Q, int d,
                                   int n) -> std::vector<uint> {
        assert((int)P.size() == d * n);
        assert((int)Q.size() == d * n);
        if (n == 1) return P;
        std::vector<uint> dftQ(d * n * 4), dftP(d * n * 4);
        for (int i = 0; i < d; ++i)
          for (int j = 0; j < n; ++j)
            dftP[(2 * (i + d) + 1) * n + j] = P[i * n + j],
                                         dftQ[i * n * 2 + j] = Q[i * n + j];
        dftQ[d * n * 2] = 1;
        FFT(d * n * 4, dftP.data(), inv_root.data());
        FFT(d * n * 4, dftQ.data(), root.data());
        P.resize(d * n * 2);
        Q.resize(d * n * 2);
        for (int i = 0; i < d * n * 4; i += 2) {
          P[i / 2] =
              ((ull)dftP[i] * dftQ[i + 1] + (ull)dftP[i + 1] * dftQ[i]) % MOD;
          if (P[i / 2] & 1) P[i / 2] += MOD;
          P[i / 2] /= 2;
          Q[i / 2] = (ull)dftQ[i] * dftQ[i + 1] % MOD;
        }
        InvFFT(d * n * 2, P.data(), root.data());
        InvFFT(d * n * 2, Q.data(), inv_root.data());
        for (int i = 0; i < d * 2; ++i)
          for (int j = 0; j < n / 2; ++j) P[i * (n / 2) + j] = P[i * n + n / 2 + j];
        assert(Q[0] == 1);
        Q[0] = 0;
        for (int i = 1; i < d * 2; ++i)
          for (int j = 0; j < n / 2; ++j) Q[i * (n / 2) + j] = Q[i * n + j];
        P.resize(d * n);
        Q.resize(d * n);
        return KinoshitaLi(KinoshitaLi, P, Q, d * 2, n / 2);
      };
      f.insert(f.begin(), len - n, 0);
      f.resize(len);
      g.resize(len);
      std::reverse(f.begin(), f.end());
      for (int i = 0; i < len; ++i) g[i] = (g[i] != 0 ? MOD - g[i] : 0);
      auto res = KinoshitaLi(KinoshitaLi, f, g, 1, len);
      res.resize(n);
      return res;
    }

    // 形式幂级数幂函数，计算 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 FPSComposition(f, g, n);
    }

    // 形式幂级数复合逆
    // 计算 g mod x^n 满足 g(f) = f(g) = x 要求 g(0) = 0 且 g'(0) ≠ 0
    std::vector<uint> FPSReversion(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 = PowerProjection(std::vector<uint>{1u}, f, n);
      for (int i = 1; i < n; ++i)
        proj[i] = (ull)proj[i] * (n - 1) % MOD * inv[i] % MOD;
      std::reverse(proj.begin(), proj.end());
      std::vector<uint> res = FPSPow1(proj, InvMod(MOD + 1 - n), n - 1);
      for (int i = 0; i < n - 1; ++i) res[i] = (ull)res[i] * invf1 % MOD;
      res.insert(res.begin(), 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 = FPSReversion(f, n);
      for (int i = 0; i < n; ++i) std::printf("%u%c", res[i], " \n"[i + 1 == n]);
      return 0;
    }

### 由转置原理导出

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

## 参考文献

1. Yasunori Kinoshita, Baitian Li.Power Series Composition in Near-Linear Time. FOCS 2024.
2. Alin Bostan, Ryuhei Mori.A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence. SOSA 2021: 118-132
3. R. P. Brent and H. T. Kung. 1978.Fast Algorithms for Manipulating Formal Power Series. J. ACM 25, 4 (Oct. 1978), 581–595.
4. 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/1/7 08:56:54，更新历史
> 发现错误？想一起完善？ 在 GitHub 上编辑此页！
> 本页面贡献者：Tiphereth-A, c-forrest, hly1204
> 本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用