浅谈分治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 存储处理较方便。代码可参考上方链接博文。