闲话
没自觉地发现活干不完
菜菜菜
昨天讲解是不是复读含量太高了?
人会约束自己。因为爱是永恒,是奉献,是克制,是善良。因此,人在发现自己越矩之时感到愧疚,更因为过去的愧疚踟蹰不前。
其实呢,其实呢,人也是可以为了私心而活的。爱也可以是冲动,是欲望,是所取。想要看到“花海盛开,燕子归来”,因为当人暂时穿越痛苦和罪恶,单纯地欣赏美的时候,事物的美丽才有它本来的归宿。
半在线卷积(分治 FFT)
突然发现我没有详细讲过这个东西的基本原理。这里提到了优化复杂度的 \(k\) 叉内容,这里则是有个实现的板子。本文主要讲一下基础的原理,照着 \(O(n\log^2 n)\) 的做法说。算是对浅谈幂级数的补充吧。
回顾多项式指数函数。已经熟知使用多项式牛顿迭代的方式求解的 \(O(n\log n)\) 的方法,其复杂度也是比较优秀。然而考虑到多项式牛顿迭代的过程中需要大量的 FFT 操作,其常数并不是很优秀。半在线卷积就提供了一种 \(O(n\log^2 n)\) 的小常数做法,在 OI 数据范围内的时间开销常常优于正常实现的牛顿迭代法。
首先我们需要了解半在线卷积所求解的对象。我们假设存在一个序列 \(\langle f_i \rangle\),其第 \(n\) 项系数按如下方式给定:
\[f_n = c_n \times \sum_{i=0}^{n-1}f_ig{n-i} \]其中序列 \(\langle g_i \rangle\) 是给定的一个序列。
我们观察每个元素 \(f_i\),可以发现其会且仅会向 \(f_j\text{ s.t. }j > i\) 的项贡献。这类特殊的性质不禁令我们想起一些关于特殊贡献方式的 DP 问题,具体地说,贡献与时间相关的一类题目。这类题目中,我们通常会抽象出多层偏序,随后用时间序列上分治的方法固定贡献对象求解。不难想到借用 CDQ 分治的思想。
我们对序列进行如 CDQ 的分治,也就是在对应的平衡树上作中序遍历。每次遍历到一个节点时,发现只有左子树对右子树的贡献。因此可以用已知的左子树和 \(\langle g_i \rangle\) 作卷积,取结果的后半部分作为左子树对右子树的贡献。
容易发现这样做的时间复杂度为 \(O(n \log^2 n)\)。
一份 P4721 的实现
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (register int i = (a), i##_ = (b) + 1; i < i##_; ++i)
#define pre(i,a,b) for (register int i = (a), i##_ = (b) - 1; i > i##_; --i)
const int N = 3e5 + 10, mod = 998244353, g = 3;
using poly = vector<int>;
int n, m;
int qp(int a, int b) {
int ret = 1;
while (b) {
if (b & 1) ret = 1ll * ret * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
} return ret;
}
const int L = 1 << 21;
int w[2][L];
int IUR() {
w[0][0] = w[1][0] = 1;
int w0 = qp(g, (mod - 1) / L);
for (int i = 1; i < L; ++ i) w[1][i] = w[0][L - i] = 1ll * w[1][i - 1] * w0 % mod;
return w0;
} int dmy = IUR();
int trs[L];
int initrs(int len) {
int limit = 1; while (limit < len) limit <<= 1;
for (int i = 0; i < limit; ++ i) trs[i] = (trs[i >> 1] >> 1) | ((i & 1) ? limit >> 1 : 0);
return limit;
}
void NTT(poly & A, int len, int type) {
A.resize(len);
for (int i = 0; i < len; ++ i) if (i < trs[i]) swap(A[i], A[trs[i]]);
for (int mid = 1; mid < len; mid <<= 1) {
for (int i = L / (mid << 1), j = 0; j < len; j += (mid << 1)) {
for (int k = 0; k < mid; ++ k) {
int x = A[j + k], y = 1ll * A[j + k + mid] * w[type][i * k] % mod;
A[j + k] = (x + y) % mod;
A[j + k + mid] = (x - y + mod) % mod;
}
}
} if (type == 1) return;
int tmp = qp(len, mod - 2);
for (int i = 0; i < len; ++ i) A[i] = 1ll * A[i] * tmp % mod;
}
void CDQ_FFT(poly & F, poly & G, int l, int r) {
if (l == r) { if (!l) F[l] = 1; return; }
int mid = l + r >> 1;
CDQ_FFT(F, G, l, mid);
static poly A, B; A.clear(), B.clear();
A.resize(mid-l+1), B.resize(r-l+1);
rep(i,l,mid) A[i-l] = F[i];
rep(i,1,r-l) B[i-1] = G[i];
int len = initrs(r - l + 1);
NTT(A, len, 1); NTT(B, len, 1);
for (int i = 0; i < len; ++ i) A[i] = 1ll * A[i] * B[i] % mod;
NTT(A, len, 0);
for (int i = mid + 1; i <= r; ++ i) F[i] = (F[i] + A[i - l - 1]) % mod;
CDQ_FFT(F, G, mid + 1, r);
}
signed main() {
cin >> n;
static poly F(n), G(n);
for (int i = 1; i < n; ++ i) cin >> G[i];
CDQ_FFT(F, G, 0, n-1);
for (int i = 0; i < n; ++ i) cout << F[i] << ' ';
}
然后我们有了这套理论就能拿去做一些模板了。依稀记得那天晚上一点多看 EI 省选课时 EI 挑幸运观众上去写 ln exp 和逆的半在线卷积形式
我们就以多项式指数函数为例。给定一个幂级数 \(F(x)\),我们知道我们需要的是 \(G(x) = \exp F(x)\)。下面就是经典的”盯着递归定义式看半个小时想出来如何比对系数“的时间。其实一般就是两边求导,有时还能得到线性递推(
对两边求导得到
\[G'(x) = F'(x) \exp F(x) = F'(x)G(x) \]这就可以比对系数了。我们对两边提取第 \(n\) 项系数得到
\[(n + 1) G[n + 1] = \sum_{i=0}^n (i + 1) F[i + 1] G[n - i] \]我们把相同项移到同一边去,得到
\[G[n] = \frac{1}{n} \sum_{i=1}^n i F[i] G[n - i] \]直接做即可。
具体的实现见 qwaszx 的博客。
标签:##_,int,闲话,mid,len,++,23.1,卷积 From: https://www.cnblogs.com/joke3579/p/chitchat230106.html