YuanAtallahRMQ1次元.cpp

#include <iostream>
#include <vector>
#include <memory>
#include <limits>
#include <random>
#include <stdint.h>
#include <functional>
#include <algorithm>

//Yuan, Hao, and Mikhail J. Atallah. "Data structures for range minimum queries in multidimensional arrays." Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms. Society for Industrial and Applied Mathematics, 2010.
//<http://yuanresearch.com/rmq_highdim.pdf>
//blockTypeは左右別々のテーブルでできるのでは？ということで、そこが論文と違うような？読み間違えてる？そもそも論文の方法だと多次元でどうなるのかわかってないのでよくわからない。
//BlockSize = 5..15 で lg blockTypes = 5,6,8,9,12,13,15,16,19,20,22
//実験してみると、candidateTableがblockTypesより全然使われてない。
//しかもBlockSizeが触れるほど使用率が減っている。
//BlockSize=14でN=10^7でも使用率≦0.5%程度。最適化の余地ありだな。TODO
template<typename Val, typename Compare = std::less<Val>, int BlockSize = 10>
class YuanAtallahRMQ {
public:
//	typedef int Val; typedef std::less<Val> Compare; static const int BlockSize = 10;

	typedef int Index; typedef char InBlockIndex; typedef unsigned BlockType;

	YuanAtallahRMQ(Compare comp_ = Compare()): comp(comp_) {
		blockHeight = 1; while((1 << blockHeight) < BlockSize) blockHeight ++;
		blockTypes = BlockType(1) << calculateMaximumTypeSize(BlockSize);

		resetBlockTables();
	}
	
	~YuanAtallahRMQ() {
		printUsedBlockTypeStatistics(minCandidateLeft.get(), "left");
		printUsedBlockTypeStatistics(minCandidateRight.get(), "right");
	}
	
	void build(const Val *a, Index n) {
		length = n;
		blocks = (n + BlockSize - 1) / BlockSize;
		treeHeight = 1; while((Index(1) << treeHeight) < blocks) ++ treeHeight;

		buildBlocksAll(a, n);
		buildBlockTreesAll(a);
	}

	//[l,r]の閉区間
	Index query(const Val *a, Index l, Index r) const {
		Index x = l / BlockSize, y = r / BlockSize, z = y - x;
		if(z == 0) return queryBlock(a, x, l % BlockSize, r % BlockSize);
		Index edge = minIndex(a,
			queryBlock(a, x, l % BlockSize, BlockSize - 1),
			queryBlock(a, y, 0, r % BlockSize));
		if(z == 1) return edge;
		if(z == 2) return minIndex(a, queryBlock(a, x+1, 0, BlockSize-1), edge);
		++ x, -- y;
		Index h = (Index)(getLCAHeight(x, y) - 1) * blocks;
		Index left = leftMin[h + y], right = rightMin[h + x];
		return minIndex(a, minIndex(a, left, right), edge);
	}
	Val queryVal(const Val *a, Index l, Index r) const {
		return a[query(a, l, r)];
	}

private:
	Compare comp;

	typedef InBlockIndex (BlockTable)[BlockSize][BlockSize];
	int blockHeight; BlockType blockTypes;
	std::unique_ptr<BlockTable[]> minCandidateLeft, minCandidateRight;

	Index length, blocks;
	int treeHeight;

	std::unique_ptr<BlockType[]> blockTypeLeft, blockTypeRight;
	std::unique_ptr<Index[]> leftMin, rightMin;
	
	template<typename T> struct RefLeft {
		const T left;
		RefLeft(T a, Index n): left(a) { }
		inline T ref(Index i) const { return left + i; }
		Index getOffset(T p) const { return p - left; }
		inline static Index queryCandidate(Index l, Index r) { return r; }
	};
	
	template<typename T> struct RefRight {
		const T right;
		RefRight(T a, Index n): right(a + n - 1) { }
		inline T ref(Index i) const { return right - i; }
		Index getOffset(T p) const { return right - p + 1; }
		inline static Index queryCandidate(Index l, Index r) { return l; }
	};

	void resetBlockTables() {
		resetBlockTable(minCandidateLeft);
		resetBlockTable(minCandidateRight);
	}

	void buildBlocksAll(const Val *a, Index n) {
		std::unique_ptr<InBlockIndex[]> tmpTree(new InBlockIndex[BlockSize * blockHeight]);
		std::unique_ptr<Val[]> tmpValues(new Val[BlockSize]);

		blockTypeLeft.reset(new BlockType[blocks]);
		buildBlocks<RefLeft>(a, n, minCandidateLeft.get(), blockTypeLeft.get(), tmpTree.get(), tmpValues.get());

		blockTypeRight.reset(new BlockType[blocks]);
		buildBlocks<RefRight>(a, n, minCandidateRight.get(), blockTypeRight.get(), tmpTree.get(), tmpValues.get());
	}

	void resetBlockTable(std::unique_ptr<BlockTable[]> &p) {
		const BlockType types = blockTypes;
		p.reset(new BlockTable[types]);
		BlockTable *q = p.get();
		for(BlockType i = 0; i < types; i ++)
			q[i][0][0] = -1;	//未初期化であるとしてマークする
	}
	
	template<template<typename T> class Ref>
	void buildBlocks(const Val *a, Index n, BlockTable blockTable[], BlockType blockTypeTable[], InBlockIndex tmpTree[], Val tmpValues[]) {
		for(Index i = 0; i < blocks; i ++) {
			const Val *block = getBlock(a, n, i, tmpValues);
			blockTypeTable[i] = buildBlock<Ref>(block, blockTable, tmpTree);
		}
	}

	//端っこのブロック用。コピコンとか要求するのがあれなので、配列に余分なBlockSizeの領域を用意しておくことで回避も可能
	const Val *getBlock(const Val *a, Index n, Index i, Val tmpValues[]) const {
		Index offset = i * BlockSize;
		if(offset + BlockSize <= n)
			return a + offset;
		for(Index j = offset; j < n; j ++)
			tmpValues[j - offset] = a[j];
		Val maxVal = a[offset];
		for(Index j = offset+1; j < n; j ++)
			if(comp(maxVal, a[j])) maxVal = a[j];
		std::fill(tmpValues + (n - offset), tmpValues + BlockSize, maxVal);
		return tmpValues;
	}

	template<template<typename T> class Ref>
	BlockType buildBlock(const Val *a, BlockTable blockTable[], InBlockIndex tmpTree[]) {
		for(int i = 0; i < BlockSize; i ++)
			tmpTree[i] = i;
		BlockType blockType = buildTree<Ref>(a, BlockSize, blockHeight, tmpTree, true);

		BlockTable &entry = blockTable[blockType];
		bool computed = entry[0][0] != -1;
//		if(entry[0][0] != -1) return blockType;

		for(int l = 0; l < BlockSize; l ++) for(int r = l; r < BlockSize; r ++) {
			InBlockIndex t = entry[l][r];
			if(l == r) {
				entry[l][r] = tmpTree[l];
			}else {
				int h = (getLCAHeight(l, r) - 1) * BlockSize;
				entry[l][r] = tmpTree[h + Ref<Index>::queryCandidate(l, r)];
			}
			if(computed) my_assert(t == entry[l][r]);
		}
		return blockType;
	}

	void buildBlockTreesAll(const Val *a) {
		leftMin.reset(new Index[blocks * treeHeight]);
		buildBlockTree<RefLeft>(a, leftMin.get());

		rightMin.reset(new Index[blocks * treeHeight]);
		buildBlockTree<RefRight>(a, rightMin.get());
	}

	template<template<typename T> class Ref>
	void buildBlockTree(const Val *a, Index mins[]) const {
		Index n = blocks;
		for(Index i = 0; i < n; i ++)
			mins[i] = queryBlock(a, i, 0, BlockSize-1);
		buildTree<Ref>(a, n, treeHeight, mins, false);
	}

	template<template<typename T> class Ref, typename LocalIndex>
	BlockType buildTree(const Val *a, Index n, int height, LocalIndex mins[], bool getType) const {
		LocalIndex *prev = mins, *cnt = mins + n;	//最初のn要素は既に設定されている
		Index comparisons = 0; BlockType blockType = 0;

		for(int i = 1; i < height; i ++) {
			Index pm = (Index)1 << (i-1), m = pm << 1;
			for(Index j = 0; j < n; j += m) {
				Index nodeSize = std::min((Index)(n - j), m);
				Ref<const Val*> rval(a + j, nodeSize);
				Ref<LocalIndex *> rcnt(cnt + j, nodeSize), rprev(prev + j, nodeSize);
				Index leftSize = rval.getOffset(a + std::min(n, (Index)(j + pm)));

				for(Index k = 0; k < leftSize; k ++)
					*rcnt.ref(k) = *rprev.ref(k);
				if(leftSize == nodeSize) continue;
				if(leftSize == 0) {
					for(Index k = 0; k < nodeSize; k ++)
						*rcnt.ref(k) = *rprev.ref(k);
					continue;
				}
				Index l = leftSize, r = nodeSize;
				LocalIndex leftMin = *rprev.ref(leftSize - 1);
				while(r - l > 0) {
					Index mid = l + (r - l) / 2;
					Index x = *rprev.ref(mid);
					if(comp(a[x], a[leftMin]) || (!comp(a[leftMin], a[x]) && x < leftMin)) {
						r = mid;
						if(getType) blockType |= BlockType(1) << comparisons;
					}else {
						l = mid+1;
					}
					comparisons ++;
				}
				for(Index k = leftSize; k < l; k ++)
					*rcnt.ref(k) = leftMin;
				for(Index k = l; k < nodeSize; k ++)
					*rcnt.ref(k) = *rprev.ref(k);
			}
			prev = cnt; cnt += n;
		}

		if(getType) blockType |= BlockType(1) << comparisons;
		return blockType;
	}

	Index queryBlock(const Val *a, Index i, InBlockIndex l, InBlockIndex r) const {
		Index base = i * BlockSize;
		Index cl = base + minCandidateLeft[blockTypeLeft[i]][l][r];
		Index cr = base + minCandidateRight[blockTypeRight[i]][l][r];
		return minIndex(a, cl, cr);
	}

	inline Index minIndex(const Val *a, Index l, Index r) const {
		return comp(a[l], a[r]) || (!comp(a[r], a[l]) && l < r) ? l : r;
	}

	static inline int getLCAHeight(Index l, Index r) {
		Index x = l ^ r;
		int k = 0;
		if(x & 0xffff0000) x >>= 16, k += 16;
		if(x & 0x0000ff00) x >>= 8,  k += 8;
		if(x & 0x000000f0) x >>= 4,  k += 4;
		if(x & 0x0000000c) x >>= 2,  k += 2;
		return k + (x >> 1) + (x != 0);
	}

	//とりあえずgreedyなシミュ
	static int calculateMaximumTypeSize(int n) {
		int h = 1; while(1 << h < n) ++ h;
		int res = 0;
		for(int i = 1; i < h; i ++) {
			Index pm = (Index)1 << (i-1), m = pm << 1;
			for(Index j = 0; j < n; j += m) {
				Index nodeSize = std::min((Index)(n - j), m);
				Index leftSize = std::min((Index)(n - j), (Index)pm);
				if(leftSize == nodeSize) continue;
				if(leftSize == 0) continue;
				int l = std::min(leftSize, nodeSize - leftSize), r = nodeSize;
				while(r - l > 0) {
					Index mid = l + (r - l) / 2;
					if(mid - l > r - (mid+1)) {
						r = mid;
					}else {
						l = mid+1;
					}
					res ++;
				}
			}
		}
		res ++;
		return res;
	}

	void printUsedBlockTypeStatistics(const BlockTable *table, const char *label) const {
		BlockType types = blockTypes;
		BlockType used = 0;
		for(BlockType i = 0; i < types; i ++)
			used += table[i][0][0] != -1;
		std::cout << label << ": used " << used << " / " << types << " (" << used * 1. / types * 100 << "%)" << std::endl;
	}

};