bitset 求解多维数点

robinyqc

2024-07-09 22:02:59

Algo. & Theory

多维数点问题是这样的：

对于两个 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)，总复杂度

由于 S\geq w，所以时间复杂度是 O(\dfrac {n^2k}w)，空间复杂度是 O(\dfrac{Sn}z)。

代码实现

#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;
}