详细揭秘：如何发明小波矩阵

以静态区间第k kk小为例。

首先假装你会归并树。归并树是啥？其实就是对归并排序的过程建树。先建立一个线段树，再自底向上归并，求出每个节点对应的区间[ l , r ] [l,r][l,r]的所有元素构成的有序序列。

归并树可以慢速二维数点。具体地，我们要查区间[ l , r ] [l,r][l,r]有多少个≤ x \le x≤x的数，可以直接找到每个线段树节点，对节点上的序列二分一下。时间复杂度是O ( log ⁡ 2 n ) \mathcal{O}(\log^2 n)O(log2n)的。要是直接拿这玩意二分求区间第k kk小，可以做到高贵的O ( log ⁡ 2 n log ⁡ V ) \mathcal{O}(\log^2 n\log V)O(log2nlogV)复杂度，这也太菜了。

究其原因是没有利用好归并树能做二维数点。不妨换维，将下标与值域互换，相当于求：值域在[ l , r ] [l,r][l,r]内的从左到右第k kk个数。这就成功把要二分的东西放到了线段树维护的维度上，那么可以把二分换成线段树二分做到O ( log ⁡ n log ⁡ V ) \mathcal{O}(\log n\log V)O(lognlogV)。

还能再给力一点吗？注意到换维过后值域是O ( n ) \mathcal{O}(n)O(n)的了，如果能O ( n ) \mathcal{O}(n)O(n)地对同一层处理出来每个数前面有几个数，那么就能再去掉一个log ⁡ n \log nlogn。于是我们开始思考怎么干这个事情。

方便起见把线段树换成01trie \text{01trie}01trie，这样可以省掉一些讨论。现在每一层的元素个数是O ( n ) \mathcal{O}(n)O(n)，值域也是O ( n ) \mathcal{O}(n)O(n)，我们当然希望能把同一层的所有数都搓到一起。于是不妨将所有左子树扔到最前面，右子树扔到最后面，然后把序列拼起来。大概是这样：

这个大概叫“蝴蝶变换”。方便起见把左子树方向称为红色，右子树方向称为蓝色。

这个树的结构就可以扔掉了，因为我只需要知道一个位置前面有几个蓝色元素，就能直接推出它在蝴蝶变换后会到哪个地方。现在在线段树二分上的过程可以写成这样：

• 从第一层开始。
• 计算[ l , r ] [l,r][l,r]之间的红色元素数量。
• 若个数≥ k \ge k≥k，则往红色子树走。否则往蓝色子树走。

我们甚至只需要知道元素个数！所以每一层做一个前缀和即可。代码大概长这样：

inlineintkth(intl,intr,intk){intres=0;l--;for(inti=29;~i;i--){intql=b[i][l],qr=b[i][r],c=(r-qr)-(l-ql);if(c>=k)l-=ql,r-=qr;elseres|=1<<i,k-=c,l=p[i]+ql,r=p[i]+qr;}returnres;}

其中b i b_ibi​是第i ii层的前缀和，p i p_ipi​是第i ii层的红色元素个数。这样我们就做到了O ( log ⁡ V ) \mathcal{O}(\log V)O(logV)查询。

还能再再再给力一点吗？时间复杂度上很难优化了，但是空间O ( n log ⁡ V ) \mathcal{O}(n\log V)O(nlogV)还是不够看。发现我们只要算1 11的个数，可以压个位，空间复杂度做到O ( n log ⁡ V w ) = O ( n ) \mathcal{O}\left(\frac{n\log V}{w}\right)=\mathcal{O}(n)O(wnlogV​)=O(n)。

然后你就吊打主席树了。恭喜你发明了小波矩阵！

structBitset{intn;vector<pair<ull,int>>b;Bitset(){}Bitset(vector<ull>a){n=(a.size()>>6)+1,b.resize(n);for(inti=0;i<a.size();i++)b[i>>6].first|=a[i]<<(i&63);for(inti=1;i<n;i++)b[i].second=b[i-1].second+__builtin_popcountll(b[i-1].first);}inlineintask(intk){returnb[k>>6].second+__builtin_popcountll(b[k>>6].first&(1ull<<(k&63))-1);}};structWavelet{intn;array<Bitset,30>b;array<int,30>p;Wavelet(vector<int>a){n=a.size();vector<int>t(n);for(inti=29;~i;i--){vector<ull>tmp(n);intx=0,y=0;for(intj=0;j<n;j++)tmp[j]=a[j]>>i&1;for(intj=0;j<n;j++)(a[j]>>i&1?t[y++]:a[x++])=a[j];b[i]=Bitset(tmp),p[i]=x;for(intj=0;j<y;j++)a[x+j]=t[j];}}inlineintkth(intl,intr,intk){intres=0;l--;for(inti=29;~i;i--){intql=b[i].ask(l),qr=b[i].ask(r),c=(r-qr)-(l-ql);if(c>=k)l-=ql,r-=qr;elseres|=1<<i,k-=c,l=p[i]+ql,r=p[i]+qr;}returnres;}};