题解：P3810 【模板】三维偏序（陌上花开）
__staring__
·
·
题解
这是一篇二进制分组套小波矩阵的题解。
静态二维偏序
三维偏序太困难了，我不会！
不妨先来考虑一个相关的问题，静态二维偏序：
有一个 m \times m 的二维平面，初始给定其上 n 个点 (x_i,y_i)，有 q 次查询每次给出 x_0,y_0，求满足 x_i \le x_0 \wedge y_i \le y_0 的点数（n,m,q 同阶）。
这题我会！用主席树即可，预处理时空复杂度 O(n\log{m})-O(n\log{m})，单次询问复杂度 O(\log{m})。
要是我不会主席树，或者主席树被卡常，或者要求 O(n) 线性空间怎么办？
我们可以使用小波矩阵 (Wavelet Matrix)，如果你不知道小波矩阵是什么，可以查看我的 专栏 获得详细讲解，这里只做简要说明。
小波矩阵是基于值域二进制分解的静态简洁数据结构，能在 O(n\log{m})-O(\frac{n\log{m}}{\omega}) 的时空复杂度预处理后，单次 O(\log{m}) 地回答区间 k 小或区间排名等查询（若查询信息带权，空间复杂度为 O(n\log{m})）。
预处理
假设当前要构建序列 a_i 的小波矩阵 b_{u,i}。
初始取二进制最高位，令 u=\log{m}。
可以发现当 $u \rarr u-1$ 时，1 个区间只会被拆成 2 个区间，且一个第 $u$ 位全 0，另一个第 $u$ 位全 1，下文分别称为 0/1 区间，通过维护 $b_{u,i}$ 的前缀和可以快速转移。
由于 $b_{u,i}=0/1$，排序是 $O(n)$ 的，并且可以用 `bitset` 压位存储，因此预处理时空复杂度 $O(n\log{m})-O(\frac{n\log{m}}{\omega})$。
### 查询区间排名
假设当前要求 $[L,R]$ 中 $\le v$ 的元素个数。
维护序列上的区间 $[l,r]$ 表示，此区间中的元素未被统计，且区间外的元素都被统计，初始令 $l=L,r=R$，以及二进制位 $u=\log{m}$。
若 $v$ 的第 $u$ 位为 1，那么所有 0 区间的元素都加入答案，将 $[l,r]$ 转移到 1 区间，否则直接将 $[l,r]$ 转移到 0 区间，令 $u \larr u-1$。
枚举结束后后统计的是 $<v$ 的元素个数，还需将答案加上 $r-l+1$。
每次转移通过 $O(1)$ 查询 $b_{u,i}$ 的前缀和可以做到 $O(1)$，因此单次查询复杂度 $O(\log{m})$，因此如果要查询 $b_{u,i}$ 的带权前缀和，只能 $O(n\log{m})$ 空间。
### 解决问题
将点按 $(x,y)$ 的关键字排序后，构建 $y_i$ 的小波矩阵，每次查询时二分出按 $(x,y)$ 关键字比较时，$(x_k,y_k) \le (x_0,y_0)$ 的最大 $k$，在小波矩阵上查询 $[1,k]$ 中 $\le y_0$ 的元素个数即为答案，总复杂度 $O(n\log{n}+(n+q)\log{m})$。
# 静态转动态
再考虑另一个相关问题，动态二维偏序：
有一个 $m \times m$ 的二维平面，有 $n$ 次操作：插入一个点 $(x,y)$，或查询有多少已插入的点中满足 $x_i \le X \wedge y_i \le Y$。
考虑如何用静态的小波矩阵解决动态问题。
### 分块
一个比较直接的想法是，设立阈值 $B$，每有 $B$ 次插入就令其为一组，构建这 $B$ 个点的小波矩阵。
查询时，枚举已构建的 $O(\frac{n}{B})$ 组小波矩阵进行查询，而对于还未成组的 $O(B)$ 次插入，直接暴力枚举检查，总复杂度为 $O(\frac{n}{B} \cdot B\log{m} + n \cdot (\frac{n}{B} \cdot \log{m} + B))$，取 $B=\sqrt{n\log{m}}$ 平衡查询，得到最优复杂度 $O(n\log{m} + n\sqrt{n\log{m}})=O(n\sqrt{n\log{m}})$。
### 自顶向下分组
考虑优化未成组暴力查询的部分，若每 $B_2$ 个插入重构为大组，未重构为大组的 $O(B_2)$ 个插入内每 $B_1$ 个插入重构为小组，未重构为小组的 $O(B_1)$ 个插入暴力查询，最终平衡复杂度得到 $O(n^{\frac{4}{3}}\log^{\frac{2}{3}}{m})$。
再取一次阈值呢？设每 $B_1$ 个插入重构为小组，每 $B_2$ 个插入重构为中组，每 $B_3$ 个插入重构为大组，平衡后复杂度为 $O(n^{\frac{5}{4}}\log^{\frac{3}{4}}{m})$。
通过计算可以发现，若每 $B_k$ 个插入重构为最大的组，一直到每 $B_1$ 个插入重构为最小的组，复杂度为 $O(kn\log{m}(\frac{n}{\log{m}})^{\frac{1}{k}})$。
### 二进制分组（自底向上合并）
但是这样的自顶向下散块分组很不实用，太麻烦了！不如转换思路，自底向上整块合并。
也即，设若干阈值 $B_1,B_2,\dots,B_k$，每 $B_1$ 个插入重构为 1 级组，每 $B_2$ 个 1 级组重构为 2 级组，…，每 $B_i$ 个 $i-1$ 级组重构为 $i$ 级组。
发现每个 $i$ 级组内含 $\prod_{u=1}^i B_u$ 个元素，因此直接令 $B_1=B_2=\dots=B_k=2$，得到 $k=\log{n}$，这样每个插入会参与 $O(\log{n})$ 次重构，每次查询最多有 $O(\log{n})$ 组。
这就是**二进制分组**，若已插入 $k$ 个元素，$u$ 级组存在当且仅当 $k$ 的第 $u$ 位为 1，如果画出合并树，容易发现与线段树本质相同。
# 二维转三维
对于三维偏序，我们将所有元素按 $(a,b,c)$ 关键字排序并去重，记 $v_i$ 表示有多少个元素和 $i$ 相同，这样 $i$ 统计到的元素都在 $[1,i-1]$。
随后将元素 $(a_i,b_i,c_i,v_i)$ 视为，查询此时二维平面上有多少点满足 $x \le b_i \wedge y \le c_i$，之后在平面上插入 $v_i$ 个点 $(b_i,c_i)$。
这样三维偏序就被处理成了动态二维偏序。
综上，将三维偏序转为动态二维偏序，使用二进制分组套小波矩阵，我们可在 $O(n\log^2{n})-O(\frac{n\log{n}}{\omega})$ 的时空复杂度预处理后，以 $O(\log^2{n})$ 的复杂度回答单次询问。
# 扩展
如果题目操作除了插入还有删除怎么办？插入 $2^k-1$ 个点后，重复插入删除，每次操作都要重构 $2^k-1$ 个点，那复杂度就假了。
对于数点这种有可减性的信息，我们可以对删除操作另维护一组负贡献的二进制分组，查询时两组做差。
若信息不具有可减性，但是删除操作是撤销上一次插入，那么可以考虑懒惰合并，即当有 3 个 $k-1$ 级组时合并其中 2 组至 $k$ 级，这样每 $O(2^k)$ 次操作才会重构一次 $k$ 级组，均摊复杂度还是 $O(\log{n})$ 的。
如果不想用二进制分组，可以直接把小波矩阵做成动态的吗？
可以，把内层的 `bitset` 换成平衡树即可。
# 实现
小波矩阵查询时，维护的区间可写成左闭右开 $(l,r]$，这样就无需对 $l$ 进行偏移。
合并时实际上不需要清空不存在的组别，可以只重构新出现的组别。
[提交记录](https://www.luogu.com.cn/record/221308799)，这里的实现魔改了专栏中的模板，可能常数较大（857ms），空间为线性（5.79mb）。
```cpp
#include <bits/stdc++.h>
using namespace std;
namespace staring
{
using LL = long long;
using ULL = unsigned long long;
#define fir first
#define sec second
#define FOR(i,a,b) for(int i = (a), i##E = (b); i <= i##E; i ++)
#define ROF(i,a,b) for(int i = (a), i##E = (b); i >= i##E; i --)
template <typename TYPE>
int gmax(TYPE &x, const TYPE& y) {return x < y ? x = y, 1 : 0;}
template <typename TYPE>
int gmin(TYPE &x, const TYPE& y) {return y < x ? x = y, 1 : 0;}
static constexpr int SIZE = 1 << 20;
static char buffin[SIZE]{}, *pin1{}, *pin2{};
static char buffout[SIZE]{}, *pout{buffout};
#define GETC() (pin1 == pin2 && (pin2 = (pin1 = buffin) + fread(buffin, 1, SIZE, stdin), pin1 == pin2)? EOF : *pin1++)
#define PUTC(c) (pout - buffout == SIZE && (fwrite(buffout, 1, SIZE, stdout), pout = buffout), (*pout++ = c))
template <typename TYPE>
void read(TYPE &x)
{
static int signf{0}, chin{0};
x = signf = 0, chin = GETC();
while(chin < '0' || chin > '9') signf |= chin == '-', chin = GETC();
while(chin >= '0' && chin <= '9') x = (x << 3) + (x << 1) + (chin ^ 48), chin = GETC();
if(signf) x = -x;
}
template <typename TYPE>
void write(TYPE x, char ch = ' ')
{
static int stack[64]{}, top{0};
!x && PUTC('0'), x < 0 && (x = -x, PUTC('-'));
while(x) stack[top++] = x % 10, x /= 10;
while(top) PUTC(stack[--top] | 48);
if(ch) PUTC(ch);
}
}using namespace staring;
constexpr int N = 1e5 + 5;
tuple <int, int, int> t[N];
pair <int, int> lshy[N];
int lshz[N], tmp0[N], tmp1[N << 1];
int res[N];
struct bitSet
{
#define ppcll __builtin_popcountll
int n, zero;
vector <pair <ULL, int>> bit;
bitSet(int _n): n(_n), zero(_n)
{
bit.resize((n >> 6) + 1, make_pair(0, 0));
}
void clear()
{
zero = n;
fill(begin(bit), end(bit), make_pair(0, 0));
}
void set(int k)
{
--zero;
bit[k >> 6].fir |= 1ull << (k & 63);
}
void init()
{
FOR(i, 1, n >> 6)
bit[i].sec = bit[i - 1].sec + ppcll(bit[i - 1].fir);
}
int rank1(int k)
{
return bit[k >> 6].sec + ppcll(bit[k >> 6].fir << (~k & 63));
}
int rank0(int k)
{
return k - rank1(k);
}
#undef ppcll
};
struct waveletMatrix
{
int len, cnt;
vector <bitSet> bit;
void init(int _len)
{
len = _len;
bit.resize(len + 1, bitSet(1 << len));
}
void build(int tot)
{
int p = tot >> len + 1 << len + 1;
int l = p + 1, r = p + (1 << len);
sort(lshy + l, lshy + r + 1);
sort(lshz + l, lshz + r + 1);
cnt = unique(lshz + l, lshz + r + 1) - lshz - l;
FOR(i, l, l + cnt - 1) tmp1[lshz[i]] = i - l;
FOR(i, 1, 1 << len) tmp0[i] = tmp1[lshy[i + p].sec];
ROF(u, __lg(cnt), 0)
{
int x = 0, y = 0;
bit[u].clear();
FOR(i, 1, 1 << len)
if(tmp0[i] >> u & 1)
tmp1[++y] = tmp0[i], bit[u].set(i);
else tmp0[++x] = tmp0[i];
bit[u].init();
memcpy(tmp0 + x + 1, tmp1 + 1, sizeof(int[y]));
}
}
int id(int u, int k, int c)
{
return !c ? bit[u].rank0(k) : bit[u].zero + bit[u].rank1(k);
}
int ranking(int tot)
{
int p = tot >> len + 1 << len + 1, L = p + 1, R = p + (1 << len);
int l = 0, r = upper_bound(lshy + L, lshy + R + 1, lshy[tot + 1]) - lshy - L;
int v = upper_bound(lshz + L, lshz + L + cnt, lshy[tot + 1].sec) - lshz - L;
int res = 0;
ROF(u, __lg(cnt), 0)
{
int c = v >> u & 1;
if(c) res += bit[u].rank0(r) - bit[u].rank0(l);
l = id(u, l, c), r = id(u, r, c);
}
return res;
}
};
void mainSolve()
{
int n, m;
read(n), read(m);
FOR(i, 1, n)
{
auto& [a, b, c] = t[i];
read(a), read(b), read(c);
}
int d = 0;
while(1 << d < n) ++d;
vector <waveletMatrix> wave(d);
FOR(i, 0, d - 1) wave[i].init(i);
sort(t + 1, t + n + 1);
FOR(i, 1, n)
{
auto [x, y, z] = t[i];
lshy[i] = make_pair(y, z);
lshz[i] = z;
}
FOR(i, 1, n)
{
int j = i;
while(j < n && lshy[j + 1] == lshy[i]) ++j;
int cur = j - i;
for(int d = i - 1; d; d &= d - 1)
{
int u = __lg(-d & d);
cur += wave[u].ranking(i - 1);
}
res[cur] += j - i + 1;
for(int d = j; d; d &= d - 1)
{
int u = __lg(-d & d);
if((j >> u) == (i - 1 >> u)) break;
wave[u].build(j);
}
i = j;
}
FOR(i, 0, n - 1)
write(res[i], '\n');
}
int main()
{
auto fileIO = [](string file)
{
string In = file + ".in";
string Out = file + ".out";
freopen(In.c_str(), "r", stdin);
freopen(Out.c_str(), "w", stdout);
};
// fileIO("");
int testCount = 1;
// read(testCount);
while(testCount--) mainSolve();
fwrite(buffout,1,pout-buffout,stdout);
return 0;
}
```