静的矩形和（クエリ先読み）

data_structure/offline_rectangle_sum.hpp · GitHub で見る

共通ページ

重み付きの点と矩形を登録し、各半開矩形に含まれる点の重みの和をまとめて求めます。

#include "lib/data_structure/offline_rectangle_sum.hpp"

OfflineRectangleSum<> rectangles;
rectangles.add_point(1, 2, 10);
rectangles.add_point(3, 4, -2);
int id = rectangles.add_query(0, 0, 3, 5);
auto answers = rectangles.solve(); // answers[id] == 10

操作	内容・返り値	計算量
OfflineRectangleSum<X=int64_t,Y=X,Sum=int64_t>(point_capacity=0,query_capacity=0)	空の構造を作り、点数・矩形数の容量を予約	\(O(1)\)
add_point(x,y,weight)	重み付きの点を追加。同じ座標も別々に数える	償却 \(O(1)\)
add_query(left,down,right,up)	[left,right) × [down,up) を追加し、0 始まりの問い合わせ番号を返す	償却 \(O(1)\)
solve() const	登録した全点を対象に、問い合わせ順の vector<Sum> を返す	\(O((n+q)\log(n+q))\)

幅または高さが 0 以下の矩形の答えは 0 です。solve は登録内容を変更せず、繰り返し呼べます。問い合わせの後に追加した点も、次の solve では全て対象になります。

座標型 X・Y にはコピーと < による厳密弱順序が必要です。デフォルト構築できない座標型にも対応します。Sum は 0 の構築・加算・減算に対応する型とし、重みと答えが表現できるように指定します。64 bit 以下の整数型では中間値のオーバーフローを避けるため、必要に応じて 128 bit 整数を使います。それ以外の型では中間の部分和も表現できる必要があります。点数・矩形数はそれぞれ INT_MAX 以下とします。

点数を \(n\)、矩形数を \(q\) として、solve は \(O((n+q)\log(n+q))\) 時間・\(O(n+q)\) 追加空間です。保持空間は \(O(n+q)\) です。比較・算術演算・コピーを定数時間とします。

依存ライブラリ

• Fenwick Sum
• 64 bit キーの安定基数ソート

検証コード

• test/offline_rectangle_sum.test.cpp

コード

#ifndef OFFLINE_RECTANGLE_SUM_HPP
#define OFFLINE_RECTANGLE_SUM_HPP

#include "../components/fenwick_sum.hpp"
#include "../components/radix_sort.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <limits>
#include <type_traits>
#include <vector>

template <class X = std::int64_t, class Y = X, class Sum = std::int64_t>
class OfflineRectangleSum {
	using Accum = std::conditional_t<std::is_integral_v<Sum> && sizeof(Sum) <= 8, __int128_t, Sum>;
	struct Point { X x; Y y; Sum weight; int rank = 0; };
	struct Query { X left, right; Y down, up; };
	struct Event { X x; int down, up, query; bool subtract; };
	std::vector<Point> points;
	std::vector<Query> queries;
	template <class A>
	static void sweep(const std::vector<Point> &ordered, const std::vector<Event> &events, int size, std::vector<Accum> &answers) {
		FenwickSum<A> tree(size);
		std::size_t p = 0;
		for (const auto &event : events) {
			while (p < ordered.size() && ordered[p].x < event.x) {
				tree.add(ordered[p].rank, A(ordered[p].weight));
				++p;
			}
			auto value = tree.template prod_as<Accum>(event.down, event.up);
			if (event.subtract) answers[event.query] -= value;
			else answers[event.query] += value;
		}
	}

public:
	explicit OfflineRectangleSum(int point_capacity = 0, int query_capacity = 0) {
		assert(point_capacity >= 0 && query_capacity >= 0);
		points.reserve(point_capacity); queries.reserve(query_capacity);
	}
	void add_point(const X &x, const Y &y, const Sum &weight) {
		assert(points.size() < std::size_t(std::numeric_limits<int>::max()));
		points.push_back({x, y, weight});
	}
	int add_query(const X &left, const Y &down, const X &right, const Y &up) {
		assert(queries.size() < std::size_t(std::numeric_limits<int>::max()));
		queries.push_back({left, right, down, up});
		return int(queries.size() - 1);
	}
	std::vector<Sum> solve() const {
		if (points.empty() || queries.empty()) return std::vector<Sum>(queries.size());
		std::vector<Accum> answers(queries.size());
		auto ordered = points;
		radix_sort_detail::sort_by_key(ordered, [](const Point &p) -> const Y & { return p.y; });
		std::vector<Y> ys;
		ys.reserve(points.size());
		for (auto &p : ordered) {
			if (ys.empty() || ys.back() < p.y) ys.push_back(p.y);
			p.rank = int(ys.size() - 1);
		}
		radix_sort_detail::sort_by_key(ordered, [](const Point &p) -> const X & { return p.x; });
		std::vector<Event> events;
		events.reserve(2 * queries.size());
		for (int i = 0; i < int(queries.size()); ++i) {
			const auto &q = queries[i];
			if (!(q.left < q.right) || !(q.down < q.up)) continue;
			int down = std::lower_bound(ys.begin(), ys.end(), q.down) - ys.begin();
			int up = std::lower_bound(ys.begin(), ys.end(), q.up) - ys.begin();
			if (down == up) continue;
			events.push_back({q.left, down, up, i, true});
			events.push_back({q.right, down, up, i, false});
		}
		radix_sort_detail::sort_by_key(events, [](const Event &e) -> const X & { return e.x; });
		bool narrow = false;
		if constexpr (std::is_integral_v<Sum> && sizeof(Sum) <= 8) {
			Accum positive = 0, negative = 0;
			for (const auto &point : points) {
				Accum value(point.weight);
				if (value > 0) positive += value;
				else negative += value;
			}
			narrow = positive <= INT64_MAX && negative >= INT64_MIN;
		}
		if (narrow) {
			if constexpr (std::is_integral_v<Sum> && sizeof(Sum) <= 8) sweep<std::int64_t>(ordered, events, int(ys.size()), answers);
		} else sweep<Accum>(ordered, events, int(ys.size()), answers);

		if constexpr (std::is_same_v<Accum, Sum>) return answers;
		else return std::vector<Sum>(answers.begin(), answers.end());
	}
};

#endif

Back to top page