闲话 23.1.6

闲话

没自觉地发现活干不完
菜菜菜

今日推歌：《Marine Bloomin’》八王子P×鹿乃

昨天讲解是不是复读含量太高了？

推荐阅读：活下去了 by EntropyIncreaser

人会约束自己。因为爱是永恒，是奉献，是克制，是善良。因此，人在发现自己越矩之时感到愧疚，更因为过去的愧疚踟蹰不前。
其实呢，其实呢，人也是可以为了私心而活的。爱也可以是冲动，是欲望，是所取。想要看到“花海盛开，燕子归来”，因为当人暂时穿越痛苦和罪恶，单纯地欣赏美的时候，事物的美丽才有它本来的归宿。

半在线卷积（分治 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 的博客。