directなRMQ的な(手抜きなのでO(n)になってないきがする)

gistfile1.cpp

template<typename Val, int BlockSize>
class DirectRMQ {
public:
  typedef int Index;	//今のところ大きくともintを仮定している(queryとか)
	typedef char InBlockIndex;
	typedef InBlockIndex (*BlockTypeRef)[BlockSize];
	
	DirectRMQ() {
		calcBallotNumbers();
		buildInnerBlockTable();
	}

	void build(const Val *a, Index n) {
		blocks = (n + BlockSize - 1) / BlockSize;
		stHeight = 0; while(1 << stHeight < blocks) ++ stHeight;
		
		blockTypes.reset(new BlockTypeRef[blocks]);
		calcBlockTypes(a, n);
		buildInnerBlockTable(a, n);
		sparseTable.reset(new Index[blocks * stHeight]);
		buildSparseTable(a);
	}

	Index query(const Val *a, Index l, Index r) const {
		Index x = l / BlockSize, y = r / BlockSize, z = y - x;
		Index k = 0, e = 1, s;
		if(z == 0) return x * BlockSize + blockTypes[x][l % BlockSize][r % BlockSize];
		Index edge = minIndex(a,
			x * BlockSize + blockTypes[x][l % BlockSize][BlockSize-1],
			y * BlockSize + blockTypes[y][0][r % BlockSize]);
		if(z == 1) return edge;
		z -= 2;
		s = ((z & 0xffff0000) != 0) << 4; z >>= s; e <<= s; k |= s;
		s = ((z & 0x0000ff00) != 0) << 3; z >>= s; e <<= s; k |= s;
		s = ((z & 0x000000f0) != 0) << 2; z >>= s; e <<= s; k |= s;
		s = ((z & 0x0000000c) != 0) << 1; z >>= s; e <<= s; k |= s;
		s = ((z & 0x00000002) != 0) << 0; z >>= s; e <<= s; k |= s;
		return minIndex(a, edge, minIndex(a,
			sparseTable[x + 1 + blocks * k],
			sparseTable[y + blocks * k - e]));
	}

private:
	int ballotNumbers[BlockSize+1][BlockSize+1];
	std::unique_ptr<InBlockIndex[][BlockSize][BlockSize]> innerBlockTable;

	Index blocks;
	int stHeight;
	std::unique_ptr<BlockTypeRef[]> blockTypes;
	std::unique_ptr<Index[]> sparseTable;

	static inline Index minIndex(const Val *a, Index x, Index y) {
		return a[x] < a[y] || (a[x] == a[y] && x < y) ? x : y;
	}

	void buildSparseTable(const Val *a) {
		Index *b = sparseTable.get();
		if(stHeight) for(Index i = 0; i < blocks; i ++)
			b[i] = i * BlockSize + blockTypes[i][0][BlockSize-1];
		for(Index t = 1; t*2 < blocks; t *= 2) {
			std::copy(b, b + blocks, b + blocks);
			b += blocks;
			for(Index i = 0; i < blocks - t; ++ i)
				b[i] = minIndex(a, b[i], b[i + t]);
		}
	}

	void buildInnerBlockTable(const Val *a, Index n) {
		for(Index i = 0; i < blocks; i ++) {
			BlockTypeRef table = blockTypes[i];
			if(table[0][0] != -1) continue;
			const Val *p = getBlock(a, n, i);
			for(InBlockIndex left = 0; left < BlockSize; left ++) {
				Val minVal = p[left];
				InBlockIndex minIndex = left;
				for(InBlockIndex right = left; right < BlockSize; right ++) {
					if(p[right] < minVal) {
						minVal = p[right];
						minIndex = right;
					}
					table[left][right] = minIndex;
				}
			}
		}
	}

	//端っこのブロック用に関数内staticなテンポラリ配列を返す
	const Val *getBlock(const Val *a, Index n, Index i) {
		Index offset = i * BlockSize;
		if(offset + BlockSize <= n)
			return a + offset;
		else {
			static Val tmp_a[BlockSize];
			std::copy(a + offset, a + n, tmp_a);
			std::fill(tmp_a + (n - offset), tmp_a + BlockSize,
				std::numeric_limits<Val>::max());
			return tmp_a;
		}
	}

	void calcBlockTypes(const Val *a, Index n) {
		Val tmp_rp[BlockSize+1];
		for(Index i = 0; i < blocks; i ++)
			blockTypes[i] = calcBlockType(getBlock(a, n, i), tmp_rp);
	}

	BlockTypeRef calcBlockType(const Val *a, Val *rp) {
		rp[0] = std::numeric_limits<Val>::min();
		int q = BlockSize, N = 0;
		for(int i = 0; i < BlockSize; i ++) {
			while(rp[q + i - BlockSize] > a[i]) {
				N += ballotNumbers[BlockSize-i-1][q];
				q --;
			}
			rp[q + i + 1 - BlockSize] = a[i];
		}
		return innerBlockTable[N];
	}

	void calcBallotNumbers() {
		for(int p = 0; p <= BlockSize; p ++) {
			for(int q = 0; q <= BlockSize; q ++) {
				if(p == 0 && q == 0)
					ballotNumbers[p][q] = 1;
				else if(p <= q)
					ballotNumbers[p][q] =
						(q ? ballotNumbers[p][q-1] : 0) +
						(p ? ballotNumbers[p-1][q] : 0);
				else
					ballotNumbers[p][q] = 0;
			}
		}
	}

	void buildInnerBlockTable() {
		int numberOfTrees = ballotNumbers[BlockSize][BlockSize];
		innerBlockTable.reset(new InBlockIndex[numberOfTrees][BlockSize][BlockSize]);
		for(int i = 0; i < numberOfTrees; i ++)
			innerBlockTable[i][0][0] = -1;	//mark as "uninitialized"
	}

};