浅谈分治FFT

置顶 | 发表于 | 评论数：

§ 1 例题一

=> luogu 4721 【模板】分治 FFT

§ 1.1 题意

给定长度为 $n-1$ 的数组 $g[1],\ g[2],\cdots,g[n-1]$，求 $f[0],\ f[1],\cdots,\ f[n-1]$，其中

边界为 $f[0]=1$。答案对 $998244353$ 取模。

$2\leq n\leq 10^5,\ \ 0\leq g[i]\lt 998244353$。

§ 1.2 分析

令 $P_{\ l,r}$ 表示多项式 $P$ 的第 $l\sim r$ 项系数，即 $P[l],\ P[l+1],\cdots,P[r]$。

观察发现，题目要求的数组 $f$ 类似于卷积，但后面的项与前面有关，不能直接 ${\rm FFT}$。

考虑分治求解，假设已知 $f_{\,l,mid}$，我们需要求出其对 $f_{\,mid+1,r}$ 的贡献。改写定义可得

可以发现左半部分对右半部分贡献为

上式左边可看作 $f_{\,l,mid} \bigotimes g_{\,0,r-l}$（其中 $\bigotimes$ 为卷积符号）。

由于模数为 $998244353$，直接用原根 $G=3$ 的 ${\rm NTT}$ 求解即可。

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

§ 1.3 代码

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
// Let f[i] = sigma{j = 1 .. i} f[i - j] g[j]
// Given g[1 .. n - 1], find f[0 .. n - 1] % 998244353
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;
const int G = 3, MOD = 998244353, MAXN = 262150;

inline void getint(int &num){
	register int ch, neg = 0;
	while(!isdigit(ch = getchar())) if(ch == '-') neg = 1;
	num = ch & 15;
	while(isdigit(ch = getchar())) num = num * 10 + (ch & 15);
	if(neg) num = -num;
}

int n, g[MAXN], f[MAXN];

inline int fastpow(int bas, int ex){
	register int res = 1; bas %= MOD;
	while(ex) {if(ex & 1) res = res * (ll)bas % MOD; bas = bas * (ll)bas % MOD, ex >>= 1;}
	return res;
}

namespace Poly{
	int rev[MAXN], lim, invlim, s, Wn[MAXN], a[MAXN], b[MAXN];
	#define add(a, b) ((a) + (b) >= MOD ? (a) + (b) - MOD : (a) + (b))
	#define minus(a, b) ((a) < (b) ? (a) - (b) + MOD : (a) - (b))
	#define mul(a, b) ((a) * ll(b) % MOD)
	
	inline void init(const int &n){
		lim = 1, s = 0; while(lim < (n << 1)) lim <<= 1, s++; invlim = fastpow(lim, MOD - 2);
		for(register int i = 1; i < lim; i++) rev[i] = rev[i >> 1] >> 1 | (i & 1) << (s - 1);
	}
	
	inline void NTT(int *a, int f){
		for(register int i = 0; i < lim; i++) if(i < rev[i]) swap(a[i], a[rev[i]]);
		for(register int i = 1; i < lim; i <<= 1){
			Wn[0] = 1, Wn[1] = fastpow(G, (MOD - 1) / (i << 1));
			for(register int j = 2; j < i; j++) Wn[j] = mul(Wn[j - 1], Wn[1]);
			for(register int j = 0; j < lim; j += i << 1){
				register int *L = a + j, *R = L + i;
				for(register int k = 0; k < i; k++){
					const int t = mul(Wn[k], R[k]);
					R[k] = minus(L[k], t), L[k] = add(L[k], t);
				}
			}
		}
		if(f == -1){
			reverse(a + 1, a + lim);
			for(register int i = 0; i < lim; i++) a[i] = mul(a[i], invlim);
		}
	}
	
	void cdq(int *f, int *g, int l, int r){
		if(l == r) return;
		const int mid = l + r >> 1;
		cdq(f, g, l, mid), init(r - l + 1);
		memset(a, 0, sizeof(int) * lim), memset(b, 0, sizeof(int) * lim);
		for(register int i = l; i <= mid; i++) a[i - l] = f[i];
		for(register int i = 0; i <= r - l; i++) b[i] = g[i];
		NTT(a, 1), NTT(b, 1);
		for(register int i = 0; i < lim; i++) a[i] = mul(a[i], b[i]);
		NTT(a, -1);
		for(register int i = mid + 1; i <= r; i++) f[i] = add(f[i], a[i - l]);
		cdq(f, g, mid + 1, r);
	}
}

int main(){
	getint(n), f[0] = 1;
	for(register int i = 1; i < n; i++) getint(g[i]);
	Poly::cdq(f, g, 0, n - 1);
	for(register int i = 0; i < n; i++) printf("%d ", f[i]);
	return puts(""), 0;
}

§ 2 例题二

=> [EZOI][1217NOI模拟赛]math(生成函数+分治FFT+高精度) 的子问题

§ 2.1 子问题

求多项式 $\prod\limits_{i=1}^n(x+i)$ 展开后各项的值。

$n\leq 10^5$。

§ 2.2 求解

分治求解，时间复杂度 $T(n)=2\ T(\frac{n}{2})+O(n\log n)=O(n\log^2n)$。

多项式用 vector 存储处理较方便。代码可参考上方链接博文。