题解：QOJ8692 Yet Another Convolution
- 37th Petrozavodsk Programming Camp, Summer 2019
- Day 4: Oleksandr Kulkov Contest 2, Tuesday, August 27, 2019
题目描述
给你两个长为 \(n\) 的数列 \(a_i, b_j\)，对所有 \(1\leq k\leq n\)，计算 \(c_k=\max\limits_{\gcd(i, j)=k}|a_i-b_j|\)。\(n\leq 10^5\)。
题解做法
首先有一个共识：如果能在 \(T(n)\) 的时间内计算 \(c_1\)，则我们能在 \(\sum_{i=1}^nT(n/i)\) 的时间内解决原问题。所以我们只需要关心如何计算 \(c_1\)。
二分答案 \(C\)。将 \(a, b\) 分别排序，假设 \(a_i\geq b_j+C\)，那么从小到大枚举 \(a_i\)，扫描加入 \(b_j\)，将问题转化为：维护一个集合 \(S\)，支持加入一个数 \(x\)，查询 \(S\) 中是否有数和 \(y\) 互质。将询问莫反：
将 \(\sum_{x\in S, d|x}1\) 记作 \(cnt[d]\)，这个显然可以枚举因数维护。这样每次插入就是枚举因数，每次查询也是枚举因数，\(1\sim n\) 所有数只会插入一次，查询一次，复杂度明显是 \(O(n\log n)\)（调和级数）。再算上外层的二分 \(O(\log V)\) 和最外层的调和级数 \(O(\log n)\)，总复杂度 \(\sum_{i=1}^nn/i\cdot \log V\log(n/i)=O(n\log^2n\log V)\)。
神秘做法
《浅谈一种互质数对与最大公约数的维护方法》杜冠成（IOI2024 集训队论文集）
仍然只考虑计算 \(c_1\)。这次不二分答案，我们直接枚举每个 \(i\)，计算与 \(i\) 互质的 \(j\) 的 \(b_j\) 的最值。但这好像有点难算，因此考虑使用论文算法。核心想法在于，我们只需要关心 \(i\) 的质因子集合，质因子集合相同的 \(i\)，计算出来的最值是一样的。
外部维护一个集合 \(S\)，初始时是 \(1\sim n\)。然后爆搜质因子集合，记录 \(i, x\) 表示我们已经搜索完了第 \(i\) 个质数，前面选出的质因子的乘积是 \(x\)。记 \(p_k\) 表示第 \(k\) 个质数，枚举 \(j>i\) 且 \(x\cdot p_j\leq n\) 的所有 \(j\)，对于每个 \(j\)，将所有 \(p_j\) 的倍数从 \(S\) 中删除，然后递归到 \(j, x\cdot p_j\)，递归完成之后撤销刚才的删除操作。递归到边界，也就是说 \(x\) 是 \(i\) 的所有质因子的乘积。记录具体选了什么质因子 \(x=q_1\times q_2\times \cdots\times q_m\)，然后再次爆搜，搜索每个质因子的次数 \(i=q_1^{c_1}\times q_2^{c_2}\times \cdots\times q_m^{c_m}\leq n\)（次数是正整数），然后统计答案。
这里要支持：从 \(S\) 中删除一个数，查询 \(\min_{x\in S}b_x\)，撤销删除操作。提前排序，使用双向链表即可维护。
根据论文，当 \(n\leq 10^6\) 时，复杂度可以视作 \(O(n\sqrt n)\)，所以最终复杂度为 \(\sum_{i=1}^n(n/i)^{1.5}=O(n^{1.5})\)。
神奇的事情
第一种做法总复杂度 \(\sum_{i=1}^nn/i\cdot \log V\log(n/i)=O(n\log^2n\log V)\)，复杂度相比于 \(c_1\) 多了一个 \(\log n\)。
第二种做法总复杂度 \(\sum_{i=1}^n(n/i)^{1.5}=O(n^{1.5})\)，复杂度和 \(c_1\) 的一样。
经过实际计算，发现确实如此，没有算错复杂度（但是由于常数原因，第一种做法更快）。为什么会有这样的差别？次数增长了，调和级数就对它无效了，很神奇。
经过 DeepSeek 老师分析，他得到如下结论：对于和式 \(S_{b,k}(n) = \sum_{i=1}^n (n/i)^b \log^k (n/i)\)，其中 \(b, k\) 为常数且 \(k\geq 0\)，渐近行为如下：
-
当 \(b > 1\)：\(S_{b,k}(n) = \Theta(n^b \log^k n)\)。
更精确地，\(S_{b,k}(n) = n^b P_k(\log n) + O(n \log^k n)\) ，其中 \(P_k\) 是次数为 \(k\) 的多项式，主项系数为 \(\zeta(b)\) （黎曼 zeta 函数）。 -
当 \(b = 1\)：\(S_{1,k}(n) = \Theta(n (\log n)^{k+1})\)。
具体地，\(S_{1,k}(n) = \dfrac{n}{k+1} (\log n)^{k+1} + O(n \log^k n)\)。 -
当 \(b < 1\)：\(S_{b,k}(n) = \Theta(n)\)。
主项常数为 \(\dfrac{\Gamma(k+1)}{(1-b)^{k+1}}\) ，即 \(S_{b,k}(n) \sim n \cdot \dfrac{\Gamma(k+1)}{(1-b)^{k+1}}\)。
例子：\(b=1, k=1\) 是第一种做法的情况；\(b=1.5,k=0\) 是第二种做法的情况；\(b=0,k=1\) 即 \(\sum_{i=1}^n\log(n/i)=O(n)\) 是一个平时可能见到的一个重要级数。
我说调和级数发散害人不浅，有没有懂的？
更神奇的事情
DeepSeek 声称：对于非负整数 \(d\)，我们观察到如下结论，其中 \(M(d)=\sum_{i\geq 0}i^d/2^i\) 是只和 \(d\) 有关的常数，\(M(0)\sim M(5)\) 的值分别为 \(2, 2, 6, 26, 150, 1082\)（增长的比较快，但仍然是常数）。
这也能收敛？这也太神奇了吧！但于此同时，DeepSeek 也称：对于非负整数 \(d\)，
这个就没有那么神秘，因为这些项都是指数级衰减的，所以总和会被最大的一项控制。总之可以把这个结论记住。
神秘做法代码
#define NDEBUG
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define debug(...) void(0)
#define endl "\n"
#endif
#ifdef ONLINE_JUDGE
#define NF
#endif
using LL = long long;
template <int N>
struct siever {int pri[N + 1], cnt;siever() : cnt(0) {fill_n(pri, N + 1, true);pri[0] = pri[1] = false;for (int i = 1; i <= N; i++) {if (pri[i]) pri[++cnt] = i;for (int j = 1; j <= cnt && i * pri[j] <= N; j++) {pri[i * pri[j]] = false;if (i % pri[j] == 0) break;}}}
};
constexpr int N = 1e5 + 10;
int n, a[N], b[N];
siever<N> sie;
int ans, c[N], rcnt;
int L[N], R[N];
vector<pair<int&, int>> und;
int ps[N], cnt, m;
void rem(int& x) { und.emplace_back(x, x); }
void rollback(int t) { while (t < (int)und.size()) und.back().first = und.back().second, und.pop_back(); }
void del(int x) {if (c[x]++ == 0) {assert(1 <= x && x <= m);debug("del(%d)\n", x);rcnt += 1;rem(L[R[x]]), rem(R[L[x]]);L[R[x]] = L[x];R[L[x]] = R[x];}
}
void query(int y) {int x = R[0];if (x == m + 1) return ;ans = max(ans, abs(b[x] - a[y]));x = L[m + 1];assert(x > 0);ans = max(ans, abs(b[x] - a[y]));
#ifdef LOCALint lst = 1e9;for (int u = L[m + 1]; u; u = L[u]) {assert(lst >= b[u]);assert(R[L[u]] == u);assert(L[R[u]] == u);assert(!c[u]);debug("%d(%d) ", u, b[u]);lst = b[u];}debug("\n");
#endif
}
void dfs2(int now, int k) {if (k > cnt) return query(now);dfs2(now, k + 1);while (now <= m / ps[k]) dfs2(now *= ps[k], k + 1);
}
void dfs(int j, int now) {dfs2(now, 1);while (sie.pri[j] <= m / now) {int p = sie.pri[j];ps[++cnt] = p;auto t = (int)und.size();debug("t = %d\n", t);for (int i = p; i <= m; i += p) del(i);dfs(j + 1, now * p);debug("rollback(%d)\n", t);rollback(t);for (int i = p; i <= m; i += p) c[i] -= 1;j += 1, --cnt;}
}
int solve(int _m) {und.clear();ans = -1e9, m = _m;vector<int> per(m);for (int i = 0; i < m; i++) per[i] = i + 1;sort(per.begin(), per.end(), [&](int i, int j) { return b[i] < b[j]; });per.insert(per.begin(), 0);per.push_back(m + 1);assert(per.front() == 0);assert(per.back() == m + 1);for (int i = 0; i + 1 < (int)per.size(); i++) R[per[i]] = per[i + 1];for (int i = 0; i + 1 < (int)per.size(); i++) L[per[i + 1]] = per[i];for (int i = 1; i <= m; i++) c[i] = 0;dfs(1, 1);return ans;
}
int main() {
#ifndef LOCALcin.tie(nullptr)->sync_with_stdio(false);
#ifndef NF
#endif
#endifund.reserve(2e6);static int va[N], vb[N];cin >> n;for (int i = 1; i <= n; i++) cin >> va[i];for (int i = 1; i <= n; i++) cin >> vb[i];for (int d = 1; d <= n; d++) {for (int i = 1; i <= n / d; i++) a[i] = va[i * d], b[i] = vb[i * d];cout << solve(n / d) << " \n"[d == n];if (d == 1) fprintf(stderr, "rcnt = %d [d = 1]\n", rcnt);}fprintf(stderr, "rcnt = %d\n", rcnt);return 0;
}