Wavelet Matrix
(DataStructure/WaveletMatrix.hpp)

View this file on GitHub
Last update: 2021-03-18 15:32:11+09:00
Include: #include "DataStructure/WaveletMatrix.hpp"
Link: https://tkmst201.github.io/Library/DataStructure/WaveletMatrix.hpp

概要

非負整数列を扱う静的データ構造です。
各要素の値が $0$ 以上 $2^\sigma$ 未満の長さ $N$ の非負整数列に対し、事前計算 $\Theta(\sigma N)$ 、空間計算量 $N + o(N)$ で様々な計算を高速に行うことができます。

WaveletMatrix(const std::vector<value_type> & v)
- 時間 $\Theta(\sigma |$v$|)$, 空間 $|$v$| + o(|$v$|)$ 配列 v で初期化
size_t size()
- $\Theta(1)$ 非負整数列の長さ $N$ を返す
T access(size_t k)
- $\Theta(\sigma)$ $A_k$ を返す
tuple<uint32_t, uint32_t, uint32_t> rank_all(size_t l, size_t r, const T & x)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する ( $x$ より小さい値の個数, $x$ と等しい値の個数, $x$ より大きい値の個数) を返す
uint32_t rank_less_than(size_t l, size_t r, const T & x)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ より小さい値の個数を返す
uint32_t rank(size_t l, size_t r, const T & x)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ と等しい値の個数を返す
uint32_t rank_more_than(size_t l, size_t r, const T & x)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ より大きい値の個数を返す
size_t select(uint32_t k, const T & x)
- $\mathcal{O}(\sigma \log{N})$ 先頭から $k$ 番目に出現する $x$ の位置 (1-indexed) を返す
T quantile(size_t l, size_t r, uint32_t k)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ の中から昇順で $k$ 番目の値を返す
size_t range_frequency(size_t l, size_t r, T val_l, T val_r)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値の個数を返す
vector<pair<T, uint32_t>> range_min_k(size_t l, size_t r, T val_l, T val_r, uint32_t k)
- 対象となる値の種類数を $D$ として $\Theta(D + \sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値のうち、同じ値はまとめて小さい値から順に $k$ 個まで取り出し、値と出現頻度の組の配列にして返す
vector<pair<T, uint32_t>> range_max_k(size_t l, size_t r, T val_l, T val_r, uint32_t k)
- 対象となる値の種類数を $D$ として $\Theta(D + \sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値のうち、同じ値はまとめて大きい値から順に $k$ 個まで取り出し、値と出現頻度の組の配列にして返す
pair<T, bool> next_value(size_t l, size_t r, T val_l, T val_r)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値の最小値を返す
pair<T, bool> prev_value(size_t l, size_t r, T val_l, T val_r)
- $\Theta(\sigma)$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値の最大値を返す
vector<pair<T, size_t>> get_rect(size_t l, size_t r, T val_l, T val_r)
- $\Theta((r - l) \sigma \log{N})$ $A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値を昇順に、値と出現位置 (0-indexed) の組の配列にして返す

コンストラクタ

以降、$\sigma :=$ BITS とします。扱う値の範囲は $0$ 以上 $2^\sigma$ 未満です。

制約

T は int, unsigned int, long long, unsigned long long
BITS は T のビット数以下

WaveletMatrix(const std::vector<value_type> & v)

配列 v で初期化します。

制約

v のすべての要素の値は非負

計算量

時間 $\Theta(\sigma |$v$|)$
空間 $|$v$| + o(|$v$|)$

メンバ関数

以下、長さ $N$ の非負整数列 $A_0, A_1, \ldots, A_{N-1}$ を対象とします。
また、SuccintBitVector の rank の計算量を $\Theta(1)$ 、select の計算量を $\mathcal{O}(\log{N})$ と仮定しています。

size_t size()

非負整数列 $A$ の長さ $N$ を返します。

計算量

$\Theta(1)$

T access(size_t k)

$A_k$ を返します。

制約

$0 \leq i < N$

計算量

$\Theta(\sigma)$

tuple<uint32_t, uint32_t, uint32_t> rank_all(size_t l, size_t r, const T & x)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ より小さい値の個数、$x$ と等しい値の個数、$x$ より大きい値の個数を順にまとめたタプルを返します。 $l = r$ のときは $(0, 0, 0)$ を返します。

制約

$0 \leq l \leq r \leq N$
$0 \leq x < 2^\sigma$

計算量

$\Theta(\sigma)$

uint32_t rank_less_than(size_t l, size_t r, const T & x)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ より小さい値の個数を返します。 $l = r$ のときは $(0, 0, 0)$ を返します。

制約

$0 \leq l \leq r \leq N$
$0 \leq x < 2^\sigma$

計算量

$\Theta(\sigma)$

uint32_t rank(size_t l, size_t r, const T & x)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ と等しい値の個数を返します。 $l = r$ のときは $(0, 0, 0)$ を返します。

制約

$0 \leq l \leq r \leq N$
$0 \leq x < 2^\sigma$

計算量

$\Theta(\sigma)$

uint32_t rank_more_than(size_t l, size_t r, const T & x)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する $x$ より大きい値の個数を返します。 $l = r$ のときは $(0, 0, 0)$ を返します。

制約

$0 \leq l \leq r \leq N$
$0 \leq x < 2^\sigma$

計算量

$\Theta(\sigma)$

size_t select(uint32_t k, const T & x)

先頭から $k$ 番目に出現する $x$ の位置 (1-indexed) を返します。 $k = 0$ または、そのような位置が存在しない場合は $N + 1$ を返します。

制約

$0 \leq k$
$0 \leq x < 2^\sigma$

計算量

$\mathcal{O}(\sigma \log{N})$

T quantile(size_t l, size_t r, uint32_t k)

$A_l, A_{l+1}, \ldots, A_{r-1}$ の中から昇順で $k$ 番目の値を返します。

制約

$0 \leq l < r \leq N$
$0 < k \leq r - l$

計算量

$\Theta(\sigma)$

size_t range_frequency(size_t l, size_t r, T val_l, T val_r)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値の個数を返します。 $l = r$ または val_l $=$ val_r のときは $0$ を返します。
実装の都合上 val_r $= ($ T の最大値 $) + 1$ を指定することはできないので、この場合は代わりに rank_all を使用してください。

制約

$0 \leq l \leq r \leq N$
$0 \leq$ val_l $\leq$ val_r $\leq 2^\sigma$

計算量

$\Theta(\sigma)$

vector<pair<T, uint32_t>> range_min_k(size_t l, size_t r, T val_l, T val_r, uint32_t k)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値のうち、同じ値はまとめて小さい値から順に $k$ 個まで取り出し、値と出現頻度の組の配列にして返します。 $l = r$ または val_l $=$ val_r のときは空の配列を返します。
実装の都合上 val_r $= ($ T の最大値 $) + 1$ を指定することはできないので注意してください。

制約

$0 \leq l \leq r \leq N$
$0 \leq$ val_l $\leq$ val_r $\leq 2^\sigma$

計算量

対象となる値の種類数を $D$ として $\Theta(D + \sigma)$

vector<pair<T, uint32_t>> range_max_k(size_t l, size_t r, T val_l, T val_r, uint32_t k)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値のうち、同じ値はまとめて大きい値から順に $k$ 個まで取り出し、値と出現頻度の組の配列にして返します。 $l = r$ または val_l $=$ val_r のときは空の配列を返します。
実装の都合上 val_r $= ($ T の最大値 $) + 1$ を指定することはできないので注意してください。

制約

$0 \leq l \leq r \leq N$
$0 \leq$ val_l $\leq$ val_r $< 2^\sigma - 1$

計算量

対象となる値の種類数を $D$ として $\Theta(D + \sigma)$

pair<T, bool> next_value(size_t l, size_t r, T val_l, T val_r)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値の最小値を返します。そのような値が存在するときは pair の $2$ 番目の値が $true$ に、存在しない場合は $false$ になります。
実装の都合上 val_r $= ($ T の最大値 $) + 1$ を指定することはできないので注意してください。

制約

$0 \leq l \leq r \leq N$
$0 \leq$ val_l $\leq$ val_r $\leq 2^\sigma$

計算量

$\Theta(\sigma)$

pair<T, bool> prev_value(size_t l, size_t r, T val_l, T val_r)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値の最大値を返します。そのような値が存在するときは pair の $2$ 番目の値が $true$ に、存在しない場合は $false$ になります。
実装の都合上 val_r $= ($ T の最大値 $) + 1$ を指定することはできないので注意してください。

制約

$0 \leq l \leq r \leq N$
$0 \leq$ val_l $\leq$ val_r $\leq 2^\sigma$

計算量

$\Theta(\sigma)$

vector<pair<T, size_t>> get_rect(size_t l, size_t r, T val_l, T val_r)

$A_l, A_{l+1}, \ldots, A_{r-1}$ に存在する val_l 以上 val_r 未満の値を昇順に、値と出現位置 (0-indexed) の組の配列にして返します。同じ値の順序は、出現位置の昇順となります。
$l = r$ または val_l $=$ val_r のときは空の配列を返します。
実装の都合上 val_r $= ($ T の最大値 $) + 1$ を指定することはできないので注意してください。

制約

$0 \leq l \leq r \leq N$
$0 \leq$ val_l $\leq$ val_r $\leq 2^\sigma$

計算量

$\Theta((r - l) \sigma \log{N})$

使用例

#include <bits/stdc++.h>
#include "DataStructure/WaveletMatrix.hpp"
using namespace std;

int main() {
	//                                          0, 1, 2, 3, 4, 5, 6, 7, 8, 9
	WaveletMatrix<int, 3> wm({3, 0, 2, 3, 4, 6, 5, 7, 1, 1});
	cout << "size = " << wm.size() << endl; // 10
	
	// 3 0 2 3 4 6 5 7 1 1
	for (int i = 0; i < wm.size(); ++i) cout << wm.access(i) << " \n"[i + 1 == wm.size()];
	
	auto [lt, r, gt] = wm.rank_all(0, 7, 3); // target = [3, 0, 2, 3, 4, 6, 5]
	cout << "< 3 (ans " << lt << "), = 3 (ans " << r << "), > 3 (ans " << gt << ")" << endl; // 2, 2, 3
	
	cout << "rank_less_than(5, 8, 7) = " << wm.rank_less_than(5, 8, 7) << endl; // 2
	cout << "rank(5, 8, 7) = " << wm.rank(5, 8, 7) << endl; // 1
	cout << "rank_more_than(5, 8, 7) = " << wm.rank_more_than(5, 8, 7) << endl; // 0
	
	cout << "select(1, 6) = " << wm.select(1, 6) << endl; // 6 (1-indexed)
	cout << "select(2, 6) = " << wm.select(2, 6) << endl; // 11 (= N + 1) (1-indexed)
	cout << "select(2, 1) = " << wm.select(2, 1) << endl; // 10
	
	cout << "quantile(2, 10, 4) = " << wm.quantile(2, 10, 4) << endl; // 3 [1, 1, 2, 3, 4, 5, 6, 7](sorted)
	
	cout << "range_frequency(1, 10, 2, 4) = " << wm.range_frequency(1, 10, 2, 4) << endl; // 2 [2, 3]
	
	puts("range_min_k");
	// target = [3, 0, 2, 3, 4, 6] -> [3, 3, 4]
	// v = 3, frequency = 2
	// v = 4, frequency = 1
	for (auto [v, fq] : wm.range_min_k(0, 6, 3, 5, 2)) cout << "v = " << v << ", frequency = " << fq << endl;
	
	puts("range_max_k");
	
	// target = [5, 7, 1, 1] -> [5, 7, 1, 1]
	// v = 7, frequency = 1
	// v = 5, frequency = 1
	for (auto [v, fq] : wm.range_max_k(6, 10, 0, 8, 2)) cout << "v = " << v << ", frequency = " << fq << endl;
	
	
	puts("range_min_k");
	// target = [1, 1] -> [1, 1]
	// v = 1, frequency = 2
	for (auto [v, fq] : wm.range_min_k(8, 10, 1, 2, 3)) cout << "v = " << v << ", frequency = " << fq << endl;
	
	// idx: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9
	//  wm: 3, 0, 2, 3, 4, 6, 5, 7, 1, 1
	
	{
		// target = [0, 2, 3] -> [2, 3]
		auto [v, b] = wm.next_value(1, 4, 2, 5);
		cout << "next_value(1, 4, 2, 5) : " << "v = " << v << ", b = " << boolalpha << b << endl; // 2, true
	}
	{
		// target = [0, 2, 3, 4, 6, 5] -> [3, 4, 5]
		auto [v, b] = wm.prev_value(1, 7, 3, 6);
		cout << "next_value(1, 7, 3, 6) : " << "v = " << v << ", b = " << boolalpha << b << endl; // 5, true
	}
	{
		// target = [2, 3] -> []
		auto [v, b] = wm.prev_value(2, 4, 0, 1);
		cout << "next_value(2, 4, 0, 1) : " << "v = " << v << ", b = " << boolalpha << b << endl; // *, false
	}
	
	// target = [2, 3, 4, 6, 5, 7, 1, 1] -> [2, 3, 1, 1]
	// v = 1, idx = 8
	// v = 1, idx = 9
	// v = 2, idx = 2
	// v = 3, idx = 3
	for (auto [v, idx] : wm.get_rect(2, 10, 1, 4)) cout << "v = " << v << ", idx = " << idx << endl;
}

参考

2020/09/03: https://www.slideshare.net/pfi/ss-15916040
2020/09/04: https://miti-7.hatenablog.com/entry/2018/04/28/152259
2020/09/07: http://algoogle.hadrori.jp/algorithm/wavelet.html

Depends on

簡潔ビットベクトル (DataStructure/SuccintBitVector.hpp)

Verified with

Code

#ifndef INCLUDE_GUARD_WAVELET_MATRIX_HPP
#define INCLUDE_GUARD_WAVELET_MATRIX_HPP

#include "DataStructure/SuccintBitVector.hpp"

#include <vector>
#include <cstdint>
#include <cassert>
#include <tuple>
#include <stack>

/**
 * @brief https://tkmst201.github.io/Library/DataStructure/WaveletMatrix.hpp
 */
template<typename T, int BITS>
struct WaveletMatrix {
	static_assert(std::is_integral<T>::value == true);
	static_assert(BITS > 0);
	static_assert(BITS <= 64);
	static_assert(BITS <= 8 * sizeof(T));
	
	using size_type = std::size_t;
	using value_type = T;
	using const_reference = const value_type &;
	using bv_type = SuccintBitVector;
	
private:
	using uint32 = std::uint32_t;
	
	size_type n;
	std::vector<bv_type> bit_vector;
	std::vector<uint32> zero;
	
public:
	explicit WaveletMatrix(const std::vector<value_type> & v) {
		build(v);
	}
	
private:
	void build(const std::vector<value_type> & v) {
		assert(!v.empty());
		n = v.size();
		for (uint32 i = 0; i < n; ++i) assert((v[i] >> BITS) == 0);
		bit_vector.assign(BITS, bv_type(n));
		zero.assign(BITS, 0);
		std::vector<value_type> ord(n), nex_ord(n);
		ord = v;
		for (int i = BITS - 1; i >= 0; --i) {
			for (uint32 j = 0; j < n; ++j) {
				if (ord[j] >> i & 1) bit_vector[i].set(j);
				else nex_ord[zero[i]++] = ord[j];
			}
			bit_vector[i].build();
			if (i == 0) break;
			uint32 p1 = zero[i];
			for (uint32 j = 0; j < n; ++j) if (ord[j] >> i & 1) nex_ord[p1++] = ord[j];
			std::swap(ord, nex_ord);
		}
	}
	
public:
	size_type size() const noexcept {
		return n;
	}
	
	value_type access(size_type k) const noexcept {
		assert(k < n);
		value_type res = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			if (bit_vector[i].access(k)) {
				res |= 1ull << i;
				k = bit_vector[i].rank1(k) + zero[i];
			}
			else k = bit_vector[i].rank0(k);
		}
		return res;
	}
	
	std::tuple<uint32, uint32, uint32> rank_all(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		if (l == r) return {0, 0, 0};
		uint32 rank_lt = 0, rank_mt = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			const uint32 l1 = bit_vector[i].rank1(l), r1 = bit_vector[i].rank1(r);
			if (x >> i & 1) {
				rank_lt += (r - l) - (r1 - l1);
				l = l1 + zero[i];
				r = r1 + zero[i];
			}
			else {
				rank_mt += r1 - l1;
				l -= l1;
				r -= r1;
			}
		}
		return {rank_lt, r - l, rank_mt};
	}
	
	uint32 rank_less_than(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		return std::get<0>(rank_all(l, r, x));
	}
	
	uint32 rank(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		return std::get<1>(rank_all(l, r, x));
	}
	
	uint32 rank_more_than(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		return std::get<2>(rank_all(l, r, x));
	}
	
	size_type select(uint32 k, const_reference x) const noexcept {
		assert((x >> BITS) == 0);
		if (k == 0 || k > size()) return size() + 1;
		uint32 pos = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			if (x >> i & 1) pos = bit_vector[i].rank1(pos) + zero[i];
			else pos = bit_vector[i].rank0(pos);
		}
		pos += k;
		if (pos > size()) return size() + 1;
		for (int i = 0; i < BITS; ++i) {
			if (x >> i & 1) {
				if (pos <= zero[i]) return size() + 1;
				pos = bit_vector[i].select1(pos - zero[i]);
			}
			else pos = bit_vector[i].select0(pos);
			if (pos > size()) return size() + 1;
		}
		return pos;
	}
	
	value_type quantile(size_type l, size_type r, uint32 k) const noexcept {
		assert(l < r);
		assert(r <= size());
		assert(0 < k);
		assert(k <= r - l);
		value_type res = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			const uint32 l1 = bit_vector[i].rank1(l), r1 = bit_vector[i].rank1(r);
			const uint32 c = (r - l) - (r1 - l1);
			if (k <= c) {
				l -= l1;
				r -= r1;
			}
			else {
				l = l1 + zero[i];
				r = r1 + zero[i];
				k -= c;
				res |= 1ull << i;
			}
		}
		return res;
	}
	
	size_type range_frequency(size_type l, size_type r, value_type val_l, value_type val_r) const noexcept {
		assert(r <= size());
		assert(((val_r - 1) >> BITS) == 0);
		if (l >= r || val_l >= val_r) return 0;
		return rank_less_than(l, r, val_r) - rank_less_than(l, r, val_l);
	}
	
private:
	struct Data {
		int i;
		uint32 l, r;
		bool f;
		value_type val;
		Data(int i, uint32 l, uint32 r, bool f, const value_type & val)
			: i(i), l(l), r(r), f(f), val(val) {}
	};
	
public:
	std::vector<std::pair<value_type, uint32>> range_min_k(
		size_type l, size_type r, value_type val_l, value_type val_r, uint32 k) const {
		assert(r <= size());
		assert(((val_r - 1) >> BITS) == 0);
		if (l >= r || val_l >= val_r || k == 0) return {};
		std::vector<std::pair<value_type, uint32>> res;
		std::stack<Data> stk;
		stk.emplace(BITS - 1, l, r, true, 0);
		while (!stk.empty()) {
			const Data dat = stk.top();
			stk.pop();
			if (dat.i == -1) {
				if (dat.val < val_r) {
					res.emplace_back(dat.val, dat.r - dat.l);
					if (res.size() == k) break;
					continue;
				}
				break;
			}
			const uint32 l1 = bit_vector[dat.i].rank1(dat.l), r1 = bit_vector[dat.i].rank1(dat.r);
			const uint32 l0 = dat.l - l1, r0 = dat.r - r1;
			const bool bit = val_l >> dat.i & 1;
			if (l1 != r1) stk.emplace(dat.i - 1, l1 + zero[dat.i], r1 + zero[dat.i], dat.f & bit, dat.val | (1ull << dat.i));
			if (!(dat.f && bit) && l0 != r0) stk.emplace(dat.i - 1, l0, r0, dat.f, dat.val);
		}
		return res;
	}
	
	std::vector<std::pair<value_type, size_type>> range_max_k(
		size_type l, size_type r, value_type val_l, value_type val_r, uint32 k) const {
		assert(r <= size());
		assert(((val_r - 1) >> BITS) == 0);
		if (l >= r || val_l >= val_r || k == 0) return {};
		--val_r;
		std::vector<std::pair<value_type, size_type>> res;
		std::stack<Data> stk;
		stk.emplace(BITS - 1, l, r, true, 0);
		while (!stk.empty()) {
			const Data dat = stk.top();
			stk.pop();
			if (dat.i == -1) {
				if (dat.val >= val_l) {
					res.emplace_back(dat.val, dat.r - dat.l);
					if (res.size() == k) break;
					continue;
				}
				break;
			}
			const uint32 l1 = bit_vector[dat.i].rank1(dat.l), r1 = bit_vector[dat.i].rank1(dat.r);
			const uint32 l0 = dat.l - l1, r0 = dat.r - r1;
			const bool bit = val_r >> dat.i & 1;
			if (l0 != r0) stk.emplace(dat.i - 1, l0, r0, dat.f & !bit, dat.val);
			if (!(dat.f && !bit) && l1 != r1) stk.emplace(dat.i - 1, l1 + zero[dat.i], r1 + zero[dat.i], dat.f, dat.val | (1ull << dat.i));
		}
		return res;
	}
	
	std::pair<value_type, bool> next_value(size_type l, size_type r, value_type val_l, value_type val_r) const {
		const auto res = range_min_k(l, r, val_l, val_r, 1);
		if (res.empty()) return {0, false};
		return {res[0].first, true};
	}
	
	std::pair<value_type, bool> prev_value(size_type l, size_type r, value_type val_l, value_type val_r) const {
		const auto res = range_max_k(l, r, val_l, val_r, 1);
		if (res.empty()) return {0, false};
		return {res[0].first, true};
	}
	
	std::vector<std::pair<value_type, size_type>> get_rect(size_type l, size_type r, value_type val_l, value_type val_r) const {
		auto points = range_min_k(l, r, val_l, val_r, r - l);
		std::vector<std::pair<value_type, size_type>> res;
		for (const auto & p : points) {
			const uint32 c = rank(0, l, p.first);
			for (uint32 i = 0; i < p.second; ++i) res.emplace_back(p.first, select(c + i + 1, p.first) - 1);
		}
		return res;
	}
};

#endif // INCLUDE_GUARD_WAVELET_MATRIX_HPP

#line 1 "DataStructure/WaveletMatrix.hpp"



#line 1 "DataStructure/SuccintBitVector.hpp"



#include <vector>
#include <cstdint>
#include <algorithm>
#include <cassert>

/**
 * @brief https://tkmst201.github.io/Library/DataStructure/SuccintBitVector.hpp
 */
struct SuccintBitVector {
	using size_type = std::size_t;
	using uint32 = std::uint32_t;
	
private:
	using uint8 = std::uint8_t;
	
	static constexpr size_type LARGE = 8;
	static constexpr size_type SMALL = 3;
	
private:
	size_type n;
	std::vector<uint8> bits;
	std::vector<uint32> large;
	std::vector<uint8> small;
	
public:
	explicit SuccintBitVector(size_type n)
		: n(n), bits(n == 0 ? 0 : ((n - 1) >> SMALL) + 1, 0u) {}
	
	explicit SuccintBitVector(const std::vector<uint8> & bits)
		: n(bits.size() << SMALL), bits(bits) {}
	
	void build() {
		large.assign(((n - 1) >> LARGE) + 1, 0);
		small.assign(((n - 1) >> SMALL) + 1, 0);
		for (uint32 i = 0, lidx = 0; i < bits.size(); i += 1u << (LARGE - SMALL), ++lidx) {
			if (lidx > 0) large[lidx] = large[lidx - 1] + small[i - 1] + pop_count(bits[i - 1]);
			small[i] = 0;
			for (uint32 j = i + 1; j < std::min<uint32>(bits.size(), i + (1u << (LARGE - SMALL))); ++j) {
				small[j] = small[j - 1] + pop_count(bits[j - 1]);
			}
		}
	}
	
	size_type size() const noexcept {
		return n;
	}
	
	void set(size_type i) noexcept {
		assert(i < n);
		bits[i >> SMALL] |= 1u << (i & ~(~0u << SMALL));
	}
	
	void reset(size_type i) noexcept {
		assert(i < n);
		bits[i >> SMALL] &= ~(1u << (i & ~(~0u << SMALL)));
	}
	
	bool access(size_type i) const noexcept {
		assert(i < n);
		return bits[i >> SMALL] >> (i & ~(~0u << SMALL)) & 1;
	}
	
	uint32 rank1(size_type i) const noexcept {
		assert(i <= n);
		if (i == 0) return 0;
		--i;
		return large[i >> LARGE] + small[i >> SMALL] + pop_count(bits[i >> SMALL] & ~(~0u << ((i & (~(~0u << SMALL))) + 1)));
	}
	
	uint32 rank0(size_type i) const noexcept {
		assert(i <= n);
		return i - rank1(i);
	}
	
	size_type select1(size_type k) const noexcept {
		if (k == 0) return size() + 1;
		uint32 l = 0, r = large.size();
		while (l + 1 < r) {
			const uint32 m = (l + r) >> 1;
			(large[m] < k ? l : r) = m;
		}
		uint32 cnt = large[l];
		l = l << (LARGE - SMALL); r = std::min<uint32>(small.size(), l + (1u << (LARGE - SMALL)));
		while (l + 1 < r) {
			const uint32 m = (l + r) >> 1;
			(cnt + small[m] < k ? l : r) = m;
		}
		cnt += small[l];
		const uint32 idx = l;
		l = 0; r = size() < l * (1u << SMALL) + 8 ? ((size() - 1) & ~(0u << SMALL)) + 1 : 8;
		if (cnt + pop_count(bits[idx] & ~(~0u << r)) < k) return size() + 1;
		while (l + 1 < r) {
			const uint32 m = (l + r) >> 1;
			(cnt + pop_count(bits[idx] & ~(~0u << m)) < k ? l : r) = m;
		}
		return (idx << SMALL) + r;
	}
	
	size_type select0(size_type k) const noexcept {
		if (k == 0) return size() + 1;
		uint32 l = 0, r = large.size();
		while (l + 1 < r) {
			const uint32 m = (l + r) >> 1;
			(m * (1 << LARGE) - large[m] < k ? l : r) = m;
		}
		uint32 cnt = l * (1 << LARGE) - large[l];
		l = l << (LARGE - SMALL); r = std::min<uint32>(small.size(), l + (1u << (LARGE - SMALL)));
		while (l + 1 < r) {
			const uint32 m = (l + r) >> 1;
			(cnt + (m & ~(~0u << (LARGE - SMALL))) * (1u << SMALL) - small[m] < k ? l : r) = m;
		}
		cnt += (l & ~(~0u << (LARGE - SMALL))) * (1u << SMALL) - small[l];
		const uint32 idx = l;
		l = 0; r = size() < l * (1u << SMALL) + 8 ? ((size() - 1) & ~(0u << SMALL)) + 1 : 8;
		if (cnt + pop_count(~bits[idx] & ~(~0u << r)) < k) return size() + 1;
		while (l + 1 < r) {
			const uint32 m = (l + r) >> 1;
			(cnt + pop_count(~bits[idx] & ~(~0u << m)) < k ? l : r) = m;
		}
		return (idx << SMALL) + r;
	}
	
private:
	static constexpr uint8 table[1u << (1u << SMALL)] {
		0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
		1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
		1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
		2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
		1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
		2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
		2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
		3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
	};
	
	static constexpr uint8 pop_count(uint8 x) noexcept {
		return table[x];
	}
};


#line 5 "DataStructure/WaveletMatrix.hpp"

#line 9 "DataStructure/WaveletMatrix.hpp"
#include <tuple>
#include <stack>

/**
 * @brief https://tkmst201.github.io/Library/DataStructure/WaveletMatrix.hpp
 */
template<typename T, int BITS>
struct WaveletMatrix {
	static_assert(std::is_integral<T>::value == true);
	static_assert(BITS > 0);
	static_assert(BITS <= 64);
	static_assert(BITS <= 8 * sizeof(T));
	
	using size_type = std::size_t;
	using value_type = T;
	using const_reference = const value_type &;
	using bv_type = SuccintBitVector;
	
private:
	using uint32 = std::uint32_t;
	
	size_type n;
	std::vector<bv_type> bit_vector;
	std::vector<uint32> zero;
	
public:
	explicit WaveletMatrix(const std::vector<value_type> & v) {
		build(v);
	}
	
private:
	void build(const std::vector<value_type> & v) {
		assert(!v.empty());
		n = v.size();
		for (uint32 i = 0; i < n; ++i) assert((v[i] >> BITS) == 0);
		bit_vector.assign(BITS, bv_type(n));
		zero.assign(BITS, 0);
		std::vector<value_type> ord(n), nex_ord(n);
		ord = v;
		for (int i = BITS - 1; i >= 0; --i) {
			for (uint32 j = 0; j < n; ++j) {
				if (ord[j] >> i & 1) bit_vector[i].set(j);
				else nex_ord[zero[i]++] = ord[j];
			}
			bit_vector[i].build();
			if (i == 0) break;
			uint32 p1 = zero[i];
			for (uint32 j = 0; j < n; ++j) if (ord[j] >> i & 1) nex_ord[p1++] = ord[j];
			std::swap(ord, nex_ord);
		}
	}
	
public:
	size_type size() const noexcept {
		return n;
	}
	
	value_type access(size_type k) const noexcept {
		assert(k < n);
		value_type res = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			if (bit_vector[i].access(k)) {
				res |= 1ull << i;
				k = bit_vector[i].rank1(k) + zero[i];
			}
			else k = bit_vector[i].rank0(k);
		}
		return res;
	}
	
	std::tuple<uint32, uint32, uint32> rank_all(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		if (l == r) return {0, 0, 0};
		uint32 rank_lt = 0, rank_mt = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			const uint32 l1 = bit_vector[i].rank1(l), r1 = bit_vector[i].rank1(r);
			if (x >> i & 1) {
				rank_lt += (r - l) - (r1 - l1);
				l = l1 + zero[i];
				r = r1 + zero[i];
			}
			else {
				rank_mt += r1 - l1;
				l -= l1;
				r -= r1;
			}
		}
		return {rank_lt, r - l, rank_mt};
	}
	
	uint32 rank_less_than(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		return std::get<0>(rank_all(l, r, x));
	}
	
	uint32 rank(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		return std::get<1>(rank_all(l, r, x));
	}
	
	uint32 rank_more_than(size_type l, size_type r, const_reference x) const noexcept {
		assert(l <= r);
		assert(r <= n);
		assert((x >> BITS) == 0);
		return std::get<2>(rank_all(l, r, x));
	}
	
	size_type select(uint32 k, const_reference x) const noexcept {
		assert((x >> BITS) == 0);
		if (k == 0 || k > size()) return size() + 1;
		uint32 pos = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			if (x >> i & 1) pos = bit_vector[i].rank1(pos) + zero[i];
			else pos = bit_vector[i].rank0(pos);
		}
		pos += k;
		if (pos > size()) return size() + 1;
		for (int i = 0; i < BITS; ++i) {
			if (x >> i & 1) {
				if (pos <= zero[i]) return size() + 1;
				pos = bit_vector[i].select1(pos - zero[i]);
			}
			else pos = bit_vector[i].select0(pos);
			if (pos > size()) return size() + 1;
		}
		return pos;
	}
	
	value_type quantile(size_type l, size_type r, uint32 k) const noexcept {
		assert(l < r);
		assert(r <= size());
		assert(0 < k);
		assert(k <= r - l);
		value_type res = 0;
		for (int i = BITS - 1; i >= 0; --i) {
			const uint32 l1 = bit_vector[i].rank1(l), r1 = bit_vector[i].rank1(r);
			const uint32 c = (r - l) - (r1 - l1);
			if (k <= c) {
				l -= l1;
				r -= r1;
			}
			else {
				l = l1 + zero[i];
				r = r1 + zero[i];
				k -= c;
				res |= 1ull << i;
			}
		}
		return res;
	}
	
	size_type range_frequency(size_type l, size_type r, value_type val_l, value_type val_r) const noexcept {
		assert(r <= size());
		assert(((val_r - 1) >> BITS) == 0);
		if (l >= r || val_l >= val_r) return 0;
		return rank_less_than(l, r, val_r) - rank_less_than(l, r, val_l);
	}
	
private:
	struct Data {
		int i;
		uint32 l, r;
		bool f;
		value_type val;
		Data(int i, uint32 l, uint32 r, bool f, const value_type & val)
			: i(i), l(l), r(r), f(f), val(val) {}
	};
	
public:
	std::vector<std::pair<value_type, uint32>> range_min_k(
		size_type l, size_type r, value_type val_l, value_type val_r, uint32 k) const {
		assert(r <= size());
		assert(((val_r - 1) >> BITS) == 0);
		if (l >= r || val_l >= val_r || k == 0) return {};
		std::vector<std::pair<value_type, uint32>> res;
		std::stack<Data> stk;
		stk.emplace(BITS - 1, l, r, true, 0);
		while (!stk.empty()) {
			const Data dat = stk.top();
			stk.pop();
			if (dat.i == -1) {
				if (dat.val < val_r) {
					res.emplace_back(dat.val, dat.r - dat.l);
					if (res.size() == k) break;
					continue;
				}
				break;
			}
			const uint32 l1 = bit_vector[dat.i].rank1(dat.l), r1 = bit_vector[dat.i].rank1(dat.r);
			const uint32 l0 = dat.l - l1, r0 = dat.r - r1;
			const bool bit = val_l >> dat.i & 1;
			if (l1 != r1) stk.emplace(dat.i - 1, l1 + zero[dat.i], r1 + zero[dat.i], dat.f & bit, dat.val | (1ull << dat.i));
			if (!(dat.f && bit) && l0 != r0) stk.emplace(dat.i - 1, l0, r0, dat.f, dat.val);
		}
		return res;
	}
	
	std::vector<std::pair<value_type, size_type>> range_max_k(
		size_type l, size_type r, value_type val_l, value_type val_r, uint32 k) const {
		assert(r <= size());
		assert(((val_r - 1) >> BITS) == 0);
		if (l >= r || val_l >= val_r || k == 0) return {};
		--val_r;
		std::vector<std::pair<value_type, size_type>> res;
		std::stack<Data> stk;
		stk.emplace(BITS - 1, l, r, true, 0);
		while (!stk.empty()) {
			const Data dat = stk.top();
			stk.pop();
			if (dat.i == -1) {
				if (dat.val >= val_l) {
					res.emplace_back(dat.val, dat.r - dat.l);
					if (res.size() == k) break;
					continue;
				}
				break;
			}
			const uint32 l1 = bit_vector[dat.i].rank1(dat.l), r1 = bit_vector[dat.i].rank1(dat.r);
			const uint32 l0 = dat.l - l1, r0 = dat.r - r1;
			const bool bit = val_r >> dat.i & 1;
			if (l0 != r0) stk.emplace(dat.i - 1, l0, r0, dat.f & !bit, dat.val);
			if (!(dat.f && !bit) && l1 != r1) stk.emplace(dat.i - 1, l1 + zero[dat.i], r1 + zero[dat.i], dat.f, dat.val | (1ull << dat.i));
		}
		return res;
	}
	
	std::pair<value_type, bool> next_value(size_type l, size_type r, value_type val_l, value_type val_r) const {
		const auto res = range_min_k(l, r, val_l, val_r, 1);
		if (res.empty()) return {0, false};
		return {res[0].first, true};
	}
	
	std::pair<value_type, bool> prev_value(size_type l, size_type r, value_type val_l, value_type val_r) const {
		const auto res = range_max_k(l, r, val_l, val_r, 1);
		if (res.empty()) return {0, false};
		return {res[0].first, true};
	}
	
	std::vector<std::pair<value_type, size_type>> get_rect(size_type l, size_type r, value_type val_l, value_type val_r) const {
		auto points = range_min_k(l, r, val_l, val_r, r - l);
		std::vector<std::pair<value_type, size_type>> res;
		for (const auto & p : points) {
			const uint32 c = rank(0, l, p.first);
			for (uint32 i = 0; i < p.second; ++i) res.emplace_back(p.first, select(c + i + 1, p.first) - 1);
		}
		return res;
	}
};

Wavelet Matrix (DataStructure/WaveletMatrix.hpp)

概要

コンストラクタ

WaveletMatrix(const std::vector<value_type> & v)

メンバ関数

size_t size()

T access(size_t k)

tuple<uint32_t, uint32_t, uint32_t> rank_all(size_t l, size_t r, const T & x)

uint32_t rank_less_than(size_t l, size_t r, const T & x)

uint32_t rank(size_t l, size_t r, const T & x)

uint32_t rank_more_than(size_t l, size_t r, const T & x)

size_t select(uint32_t k, const T & x)

T quantile(size_t l, size_t r, uint32_t k)

size_t range_frequency(size_t l, size_t r, T val_l, T val_r)

vector<pair<T, uint32_t>> range_min_k(size_t l, size_t r, T val_l, T val_r, uint32_t k)

vector<pair<T, uint32_t>> range_max_k(size_t l, size_t r, T val_l, T val_r, uint32_t k)

pair<T, bool> next_value(size_t l, size_t r, T val_l, T val_r)

pair<T, bool> prev_value(size_t l, size_t r, T val_l, T val_r)

vector<pair<T, size_t>> get_rect(size_t l, size_t r, T val_l, T val_r)

使用例

参考

Depends on

Verified with

Code

Wavelet Matrix
(DataStructure/WaveletMatrix.hpp)