bitset 求解多维数点
robinyqc
2024-07-09 22:02:59
多维数点问题是这样的：
> 对于两个 $k$ 元组 $a = (a_0, a_1, \dots, a_{k - 1})$ 与 $b = (b_0, b_1, \dots,
> b_{k - 1})$，定义 $dom(a, b) = \prod\limits_{i = 0}^{k - 1}[a_i \geq b_i]$，
> 也就是若每个 $a_i$ 都大于等于 $b_i$（$a$ 偏序 $b$） 则 $dom(a, b) = 1$，否则为 $0$。
>
> 给定 $n$ 个 $k$ 元组 $t_0, t_1, \dots, p_{n - 1}$，对于每个 $i \in [0, n)$，求 $g(i) =
> \sum\limits_{j = 0}^{n - 1} dom(t_i, t_j)$。
若 $k = 2$ 则问题就是数逆序对，直接树状数组即可。
若 $k = 3$，可以考虑 CDQ 分治，复杂度是 $O(n\log^2 n)$ 的，十分的优秀。
也可以使用 K-D 树，复杂度是 $O(n\sqrt n)$，也很不错。
若 $k = 4$，可以 CDQ 套 CDQ，或者 K-D 树，复杂度已经没那么好看了。
当 $k = 5$ 时，CDQ 实际效率显然已经不如暴力了，
K-D 树复杂度相较暴力也并未有太大优势……
不如来用 bitset 优化暴力的常数。
## 暴力做法
一个显而易见的暴力做法是直接枚举 $i, j$ 然后按每一维判断 $p_i, p_j$ 是否偏序。总复杂度 $O(n^2k)$。
一个性质是，$dom[(a_{0..i}), (b_{0..i})] = 1$
是 $dom[(a_{0..i + 1}), (b_{0..i + 1})] = 1$ 的必要条件。
那么可以交换枚举顺序，先枚举维数（元组下标）$d$，维护前 $d$ 维的偏序关系，
然后基于当前的偏序关系求解下一维。设 $b_{d, i, j}
= f(p_{i, 0..d - 1}, p_{j, 0..d - 1})$。那么：
$$
b_{d + 1, i, j} \gets b_{d, i, j} \cdot [p_{i, d} \geq p_{j, d}]
$$
实际上 $d$ 这一维空间上可以省略，因为每个 $(i, j)$ 都是独立的。也就是：
$$
b_{i, j} \gets b_{i, j} \cdot [p_{i, d} \geq p_{j, d}]
$$
总复杂度 $O(n^2k)$。
## bitset 优化
可以发现上面这个问题可以用 bitset 压位优化。开一个二维数组 `bitset<N> b[N]`。
考虑上面暴力的转移式子 $[p_{i, d} \geq p_{j, d}]$，可以想到给所有点按 $d$ 维排序，
得到 $q$，且 $q_i$ 在 $p$ 中位于 $l_i$。
从 $0$ 到 $n - 1$ 枚举 $i$，容易得到一个关系数组 $r$，$r_j$
表示 $p_{l_j} \leq q_i$。
那么只需要 $ b_{l_i} \gets \left(b_{l_i} \text{ bitand } r\right)$ 就可以了。
时间复杂度 $O(\dfrac {n^2k}w)$，空间复杂度 $O(\dfrac{n^2}{z})$。
其中 $w = 64, z = 8$，常数十分的优秀。
### 优化空间
此时空间并不是很优秀，很有可能开不下。此时可以用一个常见的手段：分组跑。
将 $S (S \geq w)$ 个点 $p_{l..l + S - 1}$ 分为一组，
每次只关心 $b_{l}\dots b_{l + S - 1}$
的值。分析复杂度:
+ 排序是 $O(nk\log n)$；
+ 求出 $r$ 每组复杂度是 $O({nk})$，
总复杂度是 $O(\dfrac{n^2k}S)$，所以上面不优化空间时求出 $r$
的复杂度是 $O(nk)$ 的；
+ 求出 $b$ 数组，每组复杂度 $O(\dfrac{Snk}w)$，总复杂度
$O(\dfrac {n^2k}w)$。
由于 $S\geq w$，所以时间复杂度是 $O(\dfrac {n^2k}w)$，
空间复杂度是 $O(\dfrac{Sn}z)$。
### 代码实现
```cpp
#include <iostream>
#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
#include <bitset>
#include <cmath>
using namespace std;
using i32 = int;
using i64 = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
template<typename T> using vec = vector<T>;
signed main()
{
u32 k, n;
cin >> k >> n;
// p[k][i].first = p[i][k], p[k][i].second = l[i]
vec<vec<pair<u32, u32>>> p(k, vec<pair<u32, u32>>(n));
for (u32 i = 0; i < n; i++) {
for (u32 d = 0; d < k; d++) {
cin >> p[d][i].first;
p[d][i].second = i;
}
}
for (u32 d = 0; d < k; d++) sort(p[d].begin(), p[d].end());
vec<u32> ans(n);
// S mentioned above.
// Note that larger bs, less time cost, but more space cost.
const u32 bs = ceil(sqrt(n)) * 16;
for (u32 i = 0; i * bs < n; i++) {
u32 l = i * bs, r = min(n, i * bs + bs);
// '80000' means n is no more than 80000.
vec<bitset<80000>> b(r - l);
for (auto &c: b) c.set();
for (u32 d = 0; d < k; d++) {
bitset<80000> s; ///< relationship array mentioned above
for (u32 j = 0, it = 0; j < n; j++) {
auto [v, q] = p[d][j];
while (it < n && p[d][it].first <= v) s.set(p[d][it++].second);
if (l <= q && q < r) b[q - l] &= s;
}
}
for (u32 j = l; j < r; j++) ans[j] = b[j - l].count();
}
for (u32 i = 0; i < n; i++) cout << ans[i] << '\n';
return 0;
}
```
相关题目：[CF1826E](https://codeforces.com/contest/1826/problem/E)。
其实在 CF 上我写了一份[英文版](https://codeforces.com/blog/entry/131374)的，欢迎来看。（求 upvote QwQ）