洛谷题库 P4721【模板】分治FFT
给定数列 g1⋯n,求解满足如下形式的数列 f:f0=1,fi=∑j=0i−1fjgi−j,∀i∈{1,2,⋯,n},对 998244353 取模。
这道题可以通过多项式求逆解决。我们在此不作讨论。
考虑其转移形式,全部都是依赖于前序已求出项的、固定系数的线性转移。则可以考虑CDQ分治。如果你愿意,也可以叫作“分治FFT”。
假设正处理区间 [l,r];令 mid=⌈(l+r)/2⌉,我们已经计算好 fl,⋯,mid。现统一处理左半区间对右半的贡献。
则令考虑的乘积项为 fi,i∈[l,mid],要向 fj,j∈(mid,r] 作贡献。有 fj←gj−ifi。故而对于 i∈[l,mid],其将分别与 gmid+1−l,⋯,r−l,gmid−l,⋯,r−l−1,⋯,g1,⋯,r−mid 作卷积。其二者下标的和即为贡献对象。故而我们构造两个多项式 F(x),G(x),有 [xi]F(x)=fi+l,[xj]G(x)=gj−1,将其二者作卷积,他们对 fj,j∈(mid,r] 的贡献即为 [xj−l−1]F(x)G(x)。
总时间复杂度 O(nlog2n)。
但我们在对 F(x),G(x) 作卷积时,因为点值与系数转化,被迫将二者调整至约 1.5(r−l+1)的长度;更别说还要向上取 2 的整次幂。结果最后我们只利用了 [xk]F(x)G(x),k∈[mid,r−mid),对整个多项式作卷积有点浪费了。同时观察发现,每一项均来源于两多项式 (mid−l+1) 个系数一一对应的乘积。有没有办法省去不必要的系数呢?
引理:有 ∑i=0n−1ωniα=n[α≡0(modn)]。
证明:当 α≡0(modn) 时显然成立:nωn0=n。
当 α≡0(modn) 时,由等比数列求和公式可得i=0∑n−1(ωnα)i=ωnα−1ωnnα−1=0。□
命题:若有多项式 F(x),G(x),deg(F(x))=n−1,deg(G(x))=m−1 (n≥m),则在仅调整 G(x) 长度至 n 的情况下对二者直接作DFT后点值相乘,将新的多项式 H(x) IDFT后,得到的系数满足
[xk]H(x)=i+j≡k (modn)∑[xi]F(x)[xj]G(x),∀k,i,j∈[0,n)∩Z
证明:记 ωnk 为 n 次单位根的 k 次幂。有 ωn=exp(2iπ/n)。令 F(x)=i=0∑n−1aixi,G(x)=j=0∑n−1bjxj,则在DFT时得到多项式F1(x)=∑i=0n−1F(ωni)xi,G1(x)=∑j=0n−1G(ωnj)xj,有
[xk]F1(x)=i=0∑n−1(ωnk)iai,[xk]G1(x)=j=0∑n−1(ωnk)jbj
则将二者点值相乘即得到
[xk]H1(x)=p=0∑n−1[xp]F1(x)[xp]G1(x)(ωnk)p=p=0∑n−1i=0∑n−1j=0∑n−1(ωnk)p(ωnp)i(ωnp)jaibj=i=0∑n−1j=0∑n−1(p=0∑n−1ωnp(k+i+j))aibj
由引理我们立刻得到
[xk]H1(x)=i=0∑n−1j=0∑n−1(n[i+j≡n−k(modn)])aibj
又因为在实现时我们对 H1(x) 作IDFT,系数从左到右反转,所以就有 [xk]H(x)=i+j≡k(modn)∑aibj。□
容易发现,由于 i,j∈[0,n),故在 modn 意义下,假如 k,i 均给定,方程 j+i≡k(modn) 的解 j 是唯一的。所以,直观来看,[xk]H(x) 是由 F(x),G(x) 的每一项循环移位相乘再相加所得到的。这就被称为循环卷积。
例如,当 n=5 时,[x2]H(x)=a0b2+a1b1+a2b0+a3b4+a4b3。
现在从这个角度,回过头来思考常规情况下FFT的本质:它也是循环卷积,我们调整其长度的原因就在于,要想获得常规意义下的卷积结果,必须让类似于上式红色部分(也即,有 i+j≥n 的解)的贡献变成 0。而在本题中,令 mid−l+1=t,则 r−l+1=2t−[2∤r−l+1]。我们只需要求解 [xk]H(x),k∈[t,2t−[2∤r−l+1])。容易发现,当 i∈[0,t) 时,有 j+i∈[0,n),故不需要担心额外的贡献。即便我们将其长度调整到 2 的次幂 2r,由于可以证明 t 一定严格小于 2r−1,故而其正确性仍然得到保证。
总而言之,我们只需要对长度为 2t−[2∤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 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 () {
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;
}
2 Responses
[…] 由卷积定理和循环卷积得到,我们对同样长为 nnn 的数列 {bn}{b_n}{bn} 作 DFTDFTDFT,并令数列 ck=DFT(a,n)kDFT(b,n)kc_k=DFT(a,n)_kDFT(b,n)_kck=DFT(a,n)kDFT(b,n)k,则有 IDFT(c,n)k=∑j+q≡k(modn)ajbq newcommandIDFT{operatorname{IDFT}}IDFT(c,n)_k=sum_{j+qequiv kpmod n}a_jb_q IDFT(c,n)k=j+q≡k(modn)∑ajbq […]
[…] 求出 F(ωnt),t=0,1,⋯ ,n−1F(w_n^t),t=0,1,cdots,n-1F(ωnt),t=0,1,⋯,n−1。则 (F(ωnt))k(F(w_n^t))^k(F(ωnt))k 就等于 F(z)F(z)F(z) 在 mod n{}bmod nmodn 意义下的 kkk 次幂求 z=ωnkz=w_n^kz=ωnk 点值的结果。这可以由单位根的性质说明。 […]