分治FFT与循环卷积

洛谷题库 P4721【模板】分治FFT

给定数列 g1ng_{1\cdots n},求解满足如下形式的数列 fff0=1,fi=j=0i1fjgij,i{1,2,,n}f_0=1,f_i=\sum_{j=0}^{i-1}f_jg_{i-j},\forall i \in \{1,2,\cdots,n\},对 998244353998244353 取模。

这道题可以通过多项式求逆解决。我们在此不作讨论。

考虑其转移形式,全部都是依赖于前序已求出项的、固定系数的线性转移。则可以考虑CDQ分治。如果你愿意,也可以叫作“分治FFT”。

假设正处理区间 [l,r][l, r];令 mid=(l+r)/2\text{mid}=\lceil (l+r)/2\rceil,我们已经计算好 fl,,midf_{l,\cdots,\text{mid}}。现统一处理左半区间对右半的贡献。

则令考虑的乘积项为 fi,i[l,mid]f_i, i \in [l, \text{mid}],要向 fj,j(mid,r]f_j, j \in (\text{mid}, r] 作贡献。有 fjgjifif_j \leftarrow g_{j-i}f_i。故而对于 i[l,mid]i \in [l, \text{mid}],其将分别与 gmid+1l,,rl,gmidl,,rl1,,g1,,rmidg_{\text{mid}+1-l, \cdots, r-l}, g_{\text{mid}-l, \cdots, r-l-1}, \cdots, g_{1, \cdots, r-\text{mid}} 作卷积。其二者下标的和即为贡献对象。故而我们构造两个多项式 F(x),G(x)F(x), G(x),有 [xi]F(x)=fi+l,[xj]G(x)=gj1[x^i]F(x)=f_{i+l}, [x^j]G(x)=g_{j-1},将其二者作卷积,他们对 fj,j(mid,r]f_j, j \in (\text{mid}, r] 的贡献即为 [xjl1]F(x)G(x)[x^{j-l-1}]F(x)G(x)

总时间复杂度 O(nlog2n)\mathrm{O}(n \log^2 n)

但我们在对 F(x),G(x)F(x), G(x) 作卷积时,因为点值与系数转化,被迫将二者调整至约 1.5(rl+1)1.5(r-l+1)的长度;更别说还要向上取 22 的整次幂。结果最后我们只利用了 [xk]F(x)G(x),k[mid,rmid)[x^k]F(x)G(x), k \in [\text{mid}, r-\text{mid}),对整个多项式作卷积有点浪费了。同时观察发现,每一项均来源于两多项式 (midl+1)(\text{mid}-l+1) 个系数一一对应的乘积。有没有办法省去不必要的系数呢?

引理:有 i=0n1ωniα=n[α0(modn)]\sum_{i=0}^{n-1}\omega_n^{i\alpha}=n[\alpha\equiv 0(\bmod n)]

证明:当 α0(modn)\alpha\equiv 0(\bmod n) 时显然成立:nωn0=nn\omega_n^0=n

α≢0(modn)\alpha \not\equiv 0(\bmod n) 时,由等比数列求和公式可得i=0n1(ωnα)i=ωnnα1ωnα1=0\sum\limits_{i=0}^{n-1}(\omega_n^\alpha)^{i}=\dfrac{\color{red}\omega_n^{n\alpha}-1}{\omega_n^\alpha-1}=\color{red}0\square

命题:若有多项式 F(x),G(x),deg(F(x))=n1,deg(G(x))=m1 (nm)F(x), G(x), \operatorname{deg}(F(x))=n-1, \operatorname{deg}(G(x))=m-1\space (n\geq m),则在仅调整 G(x)G(x) 长度至 nn 的情况下对二者直接作DFT后点值相乘,将新的多项式 H(x)H(x) IDFT后,得到的系数满足 [xk]H(x)=i+jk (modn)[xi]F(x)[xj]G(x),k,i,j[0,n)Z [x^k]H(x)=\sum_{i+j\equiv k\space (\bmod n)}[x^i]F(x)[x^j]G(x),\forall k,i,j \in [0, n)\cap \mathbb{Z}

证明:记 ωnk\omega_{n}^{k}nn 次单位根的 kk 次幂。有 ωn=exp(2iπ/n)\omega_n=\exp(2\mathrm{i}\pi/n)。令 F(x)=i=0n1aixi,G(x)=j=0n1bjxjF(x)=\sum\limits_{i=0}^{n-1}a_ix^i, G(x)=\sum\limits_{j=0}^{n-1}b_jx^j,则在DFT时得到多项式F1(x)=i=0n1F(ωni)xi,G1(x)=j=0n1G(ωnj)xjF_1(x)=\sum_{i=0}^{n-1}F(\omega_n^i)x^i,G_1(x)=\sum_{j=0}^{n-1}G(\omega_n^j)x^j,有 [xk]F1(x)=i=0n1(ωnk)iai,[xk]G1(x)=j=0n1(ωnk)jbj\left[x^k\right]F_1(x)=\sum\limits_{i=0}^{n-1}(\omega_n^k)^ia_i,\left[x^k\right]G_1(x)=\sum\limits_{j=0}^{n-1}(\omega_n^k)^jb_j 则将二者点值相乘即得到 [xk]H1(x)=p=0n1[xp]F1(x)[xp]G1(x)(ωnk)p=p=0n1i=0n1j=0n1(ωnk)p(ωnp)i(ωnp)jaibj=i=0n1j=0n1(p=0n1ωnp(k+i+j))aibj\begin{aligned} \left[x^k\right]H_1(x)& =\sum_{p=0}^{n-1}[x^p]F_1(x)[x^p]G_1(x)(\omega_n^k)^p\\ & =\sum_{p=0}^{n-1}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}(\omega_n^k)^p(\omega_n^p)^i(\omega_n^p)^ja_ib_j\\ & =\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}{\color{blue}\left(\sum_{p=0}^{n-1}\omega_n^{p(k+i+j)}\right)}a_ib_j \end{aligned}

引理我们立刻得到 [xk]H1(x)=i=0n1j=0n1(n[i+jnk(modn)])aibj\begin{aligned} \left[x^k\right]H_1(x)& =\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}{\color{blue}(n[i+j\equiv n-k(\bmod n)])}a_ib_j \end{aligned}

又因为在实现时我们对 H1(x)H_1(x) 作IDFT,系数从左到右反转,所以就有 [xk]H(x)=i+jk(modn)aibj[x^k]H(x)=\sum\limits_{i+j\equiv k(\bmod n)}a_ib_j\begin{aligned}&&\square\end{aligned}

容易发现,由于 i,j[0,n)i, j \in [0, n),故在 modn\bmod n 意义下,假如 k,ik,i 均给定,方程 j+ik(modn)j+i\equiv k(\bmod n) 的解 jj 是唯一的。所以,直观来看,[xk]H(x)[x^k]H(x) 是由 F(x),G(x)F(x), G(x) 的每一项循环移位相乘再相加所得到的。这就被称为循环卷积

例如,当 n=5n=5 时,[x2]H(x)=a0b2+a1b1+a2b0+a3b4+a4b3[x^2]H(x)=a_0b_2+a_1b_1+a_2b_0+a_3{\color{red}b_4}+a_4{\color{red}b_3}

现在从这个角度,回过头来思考常规情况下FFT的本质:它也是循环卷积,我们调整其长度的原因就在于,要想获得常规意义下的卷积结果,必须让类似于上式红色部分(也即,有 i+jni+j\geq n 的解)的贡献变成 00。而在本题中,令 midl+1=t\text{mid}-l+1=t,则 rl+1=2t[2rl+1]r-l+1=2t-[2\nmid r-l+1]。我们只需要求解 [xk]H(x),k[t,2t[2rl+1])[x^k]H(x), k\in [t,2t-[2\nmid r-l+1])。容易发现,当 i[0,t)i \in [0, t) 时,有 j+i[0,n)j+i \in [0, n),故不需要担心额外的贡献。即便我们将其长度调整到 22 的次幂 2r2^r,由于可以证明 tt 一定严格小于 2r12^{r-1},故而其正确性仍然得到保证。

总而言之,我们只需要对长度为 2t[2rl+1]2t-[2\nmid r-l+1] 的序列作DFT后点值相乘,再做IDFT,即可得到正确的系数表示。

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
  92. 92
  93. 93
  94. 94
  95. 95
  96. 96
  97. 97
#include <bits/stdc++.h> using namespace std; #define inl inline /* 快读已省略。 */ #define reint register int #define newl putchar('\n') typedef long long ll; // typedef unsigned long long ull; // typedef __int128 lll; // typedef long double llf; typedef pair <int, int> pint; #define fst first #define scd second #define all(p) begin (p), end (p) using vint = vector <int>; constexpr int N = 1<<18, P = 998244353, gen = 3, INF = 0x3f3f3f3f; int n; vint g, prep[N], f; inl ll fpow (ll a, ll b) { ll res = 1; a %= P; for (; b; b >>= 1) { if (b & 1) (res *= a) %= P; (a *= a) %= P; } return res; } inl void henkan (vint &f, int l) { static int tr[N], lst = tr[0] = 0; if (lst != l) { lst = l; for (int x = 1; x < 1<<l; ++x) tr[x] = tr[x>>1]>>1|((1<<l-1) * (x & 1)); } for (int x = 1; x < 1<<l; ++x) if (tr[x] < x) swap (f[tr[x]], f[x]); } #define clog2(x) ceil (log2 (x)) #define tomod(x) if (mod < INF) x.resize (mod, 0) inl ll inv (ll a) { return fpow (a, P - 2); } inl void NTT (vint &f, int l, bool rev) { f.resize (1<<l, 0); henkan (f, l); for (int len = 2; len <= 1<<l; len <<= 1) { const ll w_n = fpow (gen, (P - 1)/len); ll w = 1, g, h; for (int st = 0; st < 1<<l; st += len, w = 1) for (int i = 0; i < len/2; ++i, (w *= w_n) %= P) g = f[i + st], h = f[i + st + len/2] * w % P, f[i + st] = (h + g) % P, f[i + st + len/2] = (g + P - h) % P; } if (!rev) return; const ll p = inv (1<<l); for (int x = 0; x < 1<<l; ++x) f[x] = f[x] * p % P; reverse (begin (f) + 1, end (f)); } bool vis[N]; inl void proce (int len) { if (!len) return; if (vis[len]) return; vis[len] = 1; prep[len] = vint (begin (g) + 1, begin (g) + len); NTT (prep[len], clog2 (len-1), 0); proce (len>>1); proce (len+1>>1); } inl void cdq_dac (int l, int r) { if (l >= r) return; int mid = l + r >> 1, len = r - l + 1, t = clog2 (len-1); cdq_dac (l, mid); // * 使用循环卷积优化。 vint _f (begin (f) + l, begin (f) + mid + 1); NTT (_f, t, 0); for (int x = 0; x < 1<<t; ++x) _f[x] = 1ll * _f[x] * prep[len][x] % P; NTT (_f, t, 1); for (int k = mid - l; k <= r - l - 1; ++k) (f[k + l + 1] += _f[k]) %= P; cdq_dac (mid + 1, r); } int main () { /* 洛谷题库 P4721 【模板】分治 FFT 循环卷积优化 吴秋实 */ read (n); g.resize (n); f.resize (n); for (int x = 1; x < n; ++x) read (g[x]); proce (n); f[0] = 1; cdq_dac (0, n - 1); for (const int x : f) print (x), putc (' '); return 0; }
  • 2022年7月3日