rectangle_add_rectangle_sum.hpp¶
Dense grid supporting half-open rectangle additions and rectangle sum queries in O(log rows log columns), using four 2D Fenwick trees.
Verified by static_rectangle_add_rectangle_sum.
离线处理二维矩形加权与矩形求和；适合更新和查询双方都是轴对齐矩形的覆盖贡献问题。
Implementation¶
#ifndef NOYA_RECTANGLE_ADD_RECTANGLE_SUM_HPP
#define NOYA_RECTANGLE_ADD_RECTANGLE_SUM_HPP 1
/// @complexity Time: O(log R log C) per dense-grid operation, or
/// O((N + Q) log(N + Q)) for the offline sparse-coordinate algorithm.
/// Space: O(RC) for the dense structure or O(N + Q) offline.
#include "noya/fenwick_2d.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <vector>
namespace noya {
/// @brief Dense grid supporting half-open rectangle additions and rectangle
/// sum queries in O(log rows log columns), using four 2D Fenwick trees.
template <class T> struct rectangle_add_rectangle_sum {
int rows = 0;
int columns = 0;
fenwick_2d<T> constant;
fenwick_2d<T> row_weight;
fenwick_2d<T> column_weight;
fenwick_2d<T> product_weight;
rectangle_add_rectangle_sum() = default;
rectangle_add_rectangle_sum(int rows_, int columns_) {
build(rows_, columns_);
}
void build(int rows_, int columns_) {
assert(rows_ >= 0 && columns_ >= 0);
rows = rows_;
columns = columns_;
constant.build(rows + 1, columns + 1);
row_weight.build(rows + 1, columns + 1);
column_weight.build(rows + 1, columns + 1);
product_weight.build(rows + 1, columns + 1);
}
void add(int row_left, int row_right, int column_left, int column_right,
const T &delta) {
assert(0 <= row_left && row_left <= row_right && row_right <= rows);
assert(0 <= column_left && column_left <= column_right &&
column_right <= columns);
corner(row_left, column_left, delta);
corner(row_left, column_right, -delta);
corner(row_right, column_left, -delta);
corner(row_right, column_right, delta);
}
T prefix_sum(int row, int column) const {
assert(0 <= row && row <= rows);
assert(0 <= column && column <= columns);
return constant.prefix_sum(row, column) * T(row) * T(column) -
row_weight.prefix_sum(row, column) * T(column) -
column_weight.prefix_sum(row, column) * T(row) +
product_weight.prefix_sum(row, column);
}
T sum(int row_left, int row_right, int column_left, int column_right) const {
assert(0 <= row_left && row_left <= row_right && row_right <= rows);
assert(0 <= column_left && column_left <= column_right &&
column_right <= columns);
return prefix_sum(row_right, column_right) -
prefix_sum(row_left, column_right) -
prefix_sum(row_right, column_left) +
prefix_sum(row_left, column_left);
}
private:
void corner(int row, int column, const T &delta) {
constant.add(row, column, delta);
row_weight.add(row, column, delta * T(row));
column_weight.add(row, column, delta * T(column));
product_weight.add(row, column, delta * T(row) * T(column));
}
};
template <class Coordinate, class T> struct weighted_rectangle {
Coordinate left;
Coordinate down;
Coordinate right;
Coordinate up;
T weight;
};
template <class Coordinate> struct rectangle_query {
Coordinate left;
Coordinate down;
Coordinate right;
Coordinate up;
};
/// @brief Answer static rectangle-add rectangle-sum queries on sparse
/// coordinates. A rectangle is four signed corners. For a prefix ending at
/// (x,y), each dominated corner (a,b,w) contributes
/// `w (x-a) (y-b)`; expanding this product leaves four coefficient sums, all
/// maintained by one x-sweep and a Fenwick tree over compressed y-coordinates.
template <class Coordinate, class T>
std::vector<T> offline_rectangle_add_rectangle_sum(
const std::vector<weighted_rectangle<Coordinate, T>> &updates,
const std::vector<rectangle_query<Coordinate>> &queries) {
struct event {
Coordinate x;
Coordinate y;
T weight;
};
struct prefix_query {
Coordinate x;
Coordinate y;
int index;
int sign;
};
struct coefficients {
T constant{};
T x{};
T y{};
T xy{};
coefficients &operator+=(const coefficients &other) {
constant += other.constant;
x += other.x;
y += other.y;
xy += other.xy;
return *this;
}
};
std::vector<event> events;
std::vector<Coordinate> ys;
events.reserve(updates.size() * 4);
ys.reserve(updates.size() * 2);
for (const auto &update : updates) {
assert(update.left <= update.right && update.down <= update.up);
events.push_back({update.left, update.down, update.weight});
events.push_back({update.left, update.up, -update.weight});
events.push_back({update.right, update.down, -update.weight});
events.push_back({update.right, update.up, update.weight});
ys.push_back(update.down);
ys.push_back(update.up);
}
std::sort(events.begin(), events.end(),
[](const event &first, const event &second) {
return first.x < second.x;
});
std::sort(ys.begin(), ys.end());
ys.erase(std::unique(ys.begin(), ys.end()), ys.end());
std::vector<prefix_query> points;
points.reserve(queries.size() * 4);
for (int index = 0; index < int(queries.size()); index++) {
const auto &query = queries[index];
assert(query.left <= query.right && query.down <= query.up);
points.push_back({query.right, query.up, index, 1});
points.push_back({query.left, query.up, index, -1});
points.push_back({query.right, query.down, index, -1});
points.push_back({query.left, query.down, index, 1});
}
std::sort(points.begin(), points.end(),
[](const prefix_query &first, const prefix_query &second) {
return first.x < second.x;
});
std::vector<coefficients> fenwick(ys.size() + 1);
auto add = [&](int position, const coefficients &value) {
for (position++; position < int(fenwick.size());
position += position & -position) {
fenwick[position] += value;
}
};
auto prefix = [&](int right) {
coefficients result;
for (; right > 0; right -= right & -right) {
result += fenwick[right];
}
return result;
};
std::vector<T> answer(queries.size());
std::size_t next_event = 0;
for (const auto &point : points) {
while (next_event < events.size() &&
events[next_event].x <= point.x) {
const auto ¤t = events[next_event++];
T x = T(current.x);
T y = T(current.y);
add(int(std::lower_bound(ys.begin(), ys.end(), current.y) - ys.begin()),
{current.weight, current.weight * x, current.weight * y,
current.weight * x * y});
}
int dominated =
int(std::upper_bound(ys.begin(), ys.end(), point.y) - ys.begin());
coefficients sum = prefix(dominated);
T x = T(point.x);
T y = T(point.y);
T value = x * y * sum.constant - y * sum.x - x * sum.y + sum.xy;
if (point.sign > 0) {
answer[point.index] += value;
} else {
answer[point.index] -= value;
}
}
return answer;
}
} // namespace noya
#endif // NOYA_RECTANGLE_ADD_RECTANGLE_SUM_HPP
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <vector>
/// @complexity Time: O(log R log C) per dense-grid operation, or
/// O((N + Q) log(N + Q)) for the offline sparse-coordinate algorithm.
/// Space: O(RC) for the dense structure or O(N + Q) offline.
/// @complexity Time: O(log R log C) point update or rectangle query.
/// Space: O(RC).
namespace noya {
/// @brief Dense two-dimensional Fenwick tree for point additions and rectangle
/// sums.
template <class T> struct fenwick_2d {
int rows = 0;
int columns = 0;
std::vector<std::vector<T>> data;
fenwick_2d() = default;
fenwick_2d(int rows_, int columns_) { build(rows_, columns_); }
/// @brief Reset to a rows by columns grid of zeros.
void build(int rows_, int columns_) {
assert(rows_ >= 0 && columns_ >= 0);
rows = rows_;
columns = columns_;
data.assign(rows + 1, std::vector<T>(columns + 1));
}
/// @brief Add delta to cell (row, column).
void add(int row, int column, const T &delta) {
assert(0 <= row && row < rows);
assert(0 <= column && column < columns);
for (int x = row + 1; x <= rows; x += x & -x) {
for (int y = column + 1; y <= columns; y += y & -y) {
data[x][y] += delta;
}
}
}
/// @brief Return the sum over [0, row) x [0, column).
T prefix_sum(int row, int column) const {
assert(0 <= row && row <= rows);
assert(0 <= column && column <= columns);
T result{};
for (int x = row; x > 0; x -= x & -x) {
for (int y = column; y > 0; y -= y & -y) {
result += data[x][y];
}
}
return result;
}
/// @brief Return the sum over [row_left, row_right) x
/// [column_left, column_right).
T rectangle_sum(int row_left, int row_right, int column_left,
int column_right) const {
assert(0 <= row_left && row_left <= row_right && row_right <= rows);
assert(0 <= column_left && column_left <= column_right &&
column_right <= columns);
return prefix_sum(row_right, column_right) -
prefix_sum(row_left, column_right) -
prefix_sum(row_right, column_left) +
prefix_sum(row_left, column_left);
}
};
} // namespace noya
namespace noya {
/// @brief Dense grid supporting half-open rectangle additions and rectangle
/// sum queries in O(log rows log columns), using four 2D Fenwick trees.
template <class T> struct rectangle_add_rectangle_sum {
int rows = 0;
int columns = 0;
fenwick_2d<T> constant;
fenwick_2d<T> row_weight;
fenwick_2d<T> column_weight;
fenwick_2d<T> product_weight;
rectangle_add_rectangle_sum() = default;
rectangle_add_rectangle_sum(int rows_, int columns_) {
build(rows_, columns_);
}
void build(int rows_, int columns_) {
assert(rows_ >= 0 && columns_ >= 0);
rows = rows_;
columns = columns_;
constant.build(rows + 1, columns + 1);
row_weight.build(rows + 1, columns + 1);
column_weight.build(rows + 1, columns + 1);
product_weight.build(rows + 1, columns + 1);
}
void add(int row_left, int row_right, int column_left, int column_right,
const T &delta) {
assert(0 <= row_left && row_left <= row_right && row_right <= rows);
assert(0 <= column_left && column_left <= column_right &&
column_right <= columns);
corner(row_left, column_left, delta);
corner(row_left, column_right, -delta);
corner(row_right, column_left, -delta);
corner(row_right, column_right, delta);
}
T prefix_sum(int row, int column) const {
assert(0 <= row && row <= rows);
assert(0 <= column && column <= columns);
return constant.prefix_sum(row, column) * T(row) * T(column) -
row_weight.prefix_sum(row, column) * T(column) -
column_weight.prefix_sum(row, column) * T(row) +
product_weight.prefix_sum(row, column);
}
T sum(int row_left, int row_right, int column_left, int column_right) const {
assert(0 <= row_left && row_left <= row_right && row_right <= rows);
assert(0 <= column_left && column_left <= column_right &&
column_right <= columns);
return prefix_sum(row_right, column_right) -
prefix_sum(row_left, column_right) -
prefix_sum(row_right, column_left) +
prefix_sum(row_left, column_left);
}
private:
void corner(int row, int column, const T &delta) {
constant.add(row, column, delta);
row_weight.add(row, column, delta * T(row));
column_weight.add(row, column, delta * T(column));
product_weight.add(row, column, delta * T(row) * T(column));
}
};
template <class Coordinate, class T> struct weighted_rectangle {
Coordinate left;
Coordinate down;
Coordinate right;
Coordinate up;
T weight;
};
template <class Coordinate> struct rectangle_query {
Coordinate left;
Coordinate down;
Coordinate right;
Coordinate up;
};
/// @brief Answer static rectangle-add rectangle-sum queries on sparse
/// coordinates. A rectangle is four signed corners. For a prefix ending at
/// (x,y), each dominated corner (a,b,w) contributes
/// `w (x-a) (y-b)`; expanding this product leaves four coefficient sums, all
/// maintained by one x-sweep and a Fenwick tree over compressed y-coordinates.
template <class Coordinate, class T>
std::vector<T> offline_rectangle_add_rectangle_sum(
const std::vector<weighted_rectangle<Coordinate, T>> &updates,
const std::vector<rectangle_query<Coordinate>> &queries) {
struct event {
Coordinate x;
Coordinate y;
T weight;
};
struct prefix_query {
Coordinate x;
Coordinate y;
int index;
int sign;
};
struct coefficients {
T constant{};
T x{};
T y{};
T xy{};
coefficients &operator+=(const coefficients &other) {
constant += other.constant;
x += other.x;
y += other.y;
xy += other.xy;
return *this;
}
};
std::vector<event> events;
std::vector<Coordinate> ys;
events.reserve(updates.size() * 4);
ys.reserve(updates.size() * 2);
for (const auto &update : updates) {
assert(update.left <= update.right && update.down <= update.up);
events.push_back({update.left, update.down, update.weight});
events.push_back({update.left, update.up, -update.weight});
events.push_back({update.right, update.down, -update.weight});
events.push_back({update.right, update.up, update.weight});
ys.push_back(update.down);
ys.push_back(update.up);
}
std::sort(events.begin(), events.end(),
[](const event &first, const event &second) {
return first.x < second.x;
});
std::sort(ys.begin(), ys.end());
ys.erase(std::unique(ys.begin(), ys.end()), ys.end());
std::vector<prefix_query> points;
points.reserve(queries.size() * 4);
for (int index = 0; index < int(queries.size()); index++) {
const auto &query = queries[index];
assert(query.left <= query.right && query.down <= query.up);
points.push_back({query.right, query.up, index, 1});
points.push_back({query.left, query.up, index, -1});
points.push_back({query.right, query.down, index, -1});
points.push_back({query.left, query.down, index, 1});
}
std::sort(points.begin(), points.end(),
[](const prefix_query &first, const prefix_query &second) {
return first.x < second.x;
});
std::vector<coefficients> fenwick(ys.size() + 1);
auto add = [&](int position, const coefficients &value) {
for (position++; position < int(fenwick.size());
position += position & -position) {
fenwick[position] += value;
}
};
auto prefix = [&](int right) {
coefficients result;
for (; right > 0; right -= right & -right) {
result += fenwick[right];
}
return result;
};
std::vector<T> answer(queries.size());
std::size_t next_event = 0;
for (const auto &point : points) {
while (next_event < events.size() &&
events[next_event].x <= point.x) {
const auto ¤t = events[next_event++];
T x = T(current.x);
T y = T(current.y);
add(int(std::lower_bound(ys.begin(), ys.end(), current.y) - ys.begin()),
{current.weight, current.weight * x, current.weight * y,
current.weight * x * y});
}
int dominated =
int(std::upper_bound(ys.begin(), ys.end(), point.y) - ys.begin());
coefficients sum = prefix(dominated);
T x = T(point.x);
T y = T(point.y);
T value = x * y * sum.constant - y * sum.x - x * sum.y + sum.xy;
if (point.sign > 0) {
answer[point.index] += value;
} else {
answer[point.index] -= value;
}
}
return answer;
}
} // namespace noya