分治FFT与循环卷积

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

给定数列 $g_{1\cdots n}$，求解满足如下形式的数列 $f$：$f_0=1,f_i=\sum_{j=0}^{i-1}f_jg_{i-j},\forall i \in \{1,2,\cdots,n\}$，对 $998244353$ 取模。

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

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

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

则令考虑的乘积项为 $f_i, i \in [l, \text{mid}]$，要向 $f_j, j \in (\text{mid}, r]$ 作贡献。有 $f_j \leftarrow g_{j-i}f_i$。故而对于 $i \in [l, \text{mid}]$，其将分别与 $g_{\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)$，有 $[x^i]F(x)=f_{i+l}, [x^j]G(x)=g_{j-1}$，将其二者作卷积，他们对 $f_j, j \in (\text{mid}, r]$ 的贡献即为 $[x^{j-l-1}]F(x)G(x)$。

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

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

引理：有 $\sum_{i=0}^{n-1}\omega_n^{i\alpha}=n[\alpha\equiv 0(\bmod n)]$。

证明：当 $\alpha\equiv 0(\bmod n)$ 时显然成立：$n\omega_n^0=n$。

当 $\alpha \not\equiv 0(\bmod n)$ 时，由等比数列求和公式可得$\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), \operatorname{deg}(F(x))=n-1, \operatorname{deg}(G(x))=m-1\space (n\geq m)$，则在仅调整 $G(x)$ 长度至 $n$ 的情况下对二者直接作DFT后点值相乘，将新的多项式 $H(x)$ IDFT后，得到的系数满足
$$
[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}
$$

证明：记 $\omega_{n}^{k}$ 为 $n$ 次单位根的 $k$ 次幂。有 $\omega_n=\exp(2\mathrm{i}\pi/n)$。令 $F(x)=\sum\limits_{i=0}^{n-1}a_ix^i, G(x)=\sum\limits_{j=0}^{n-1}b_jx^j$，则在DFT时得到多项式$F_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$，有
$$\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$$
则将二者点值相乘即得到
$$\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}$$

由引理我们立刻得到
$$\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}$$

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

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

例如，当 $n=5$ 时，$[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+j\geq n$ 的解）的贡献变成 $0$。而在本题中，令 $\text{mid}-l+1=t$，则 $r-l+1=2t-[2\nmid r-l+1]$。我们只需要求解 $[x^k]H(x), k\in [t,2t-[2\nmid r-l+1])$。容易发现，当 $i \in [0, t)$ 时，有 $j+i \in [0, n)$，故不需要担心额外的贡献。即便我们将其长度调整到 $2$ 的次幂 $2^r$，由于可以证明 $t$ 一定严格小于 $2^{r-1}$，故而其正确性仍然得到保证。

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

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
#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;
}