Привет, codeforces!
Сегодня, ориентируясь на статью Триангуляция полигонов с сайта neerc.ifmo.ru, я реализовал триангуляцию любых многоугольников без самопересечений (в том числе не являющихся выпуклыми) на C++.
$$$\text{ }$$$ $$$\text{ }$$$
Зачем мне это понадобилось? На отбор на RuCode 7.0 в этом году дали задачу H. Территориальное размежевание. Задача ставится следующим образом: вам дана бесконечная плоскость, покрашенная зеброй (чёрная полоска высоты $$$1$$$, белая полоска высоты $$$1$$$, чёрная ..., белая ...). На плоскости дан многоугольник без самопересечений, состоящий из $$$n\text{ }\left(3\le n \le 100000\right)$$$ вершин с целочисленными координатами $$$x_i, y_i \left(-10^9 \le x_i, y_i \le 10^9\right)$$$. Вершины указаны в порядке какого-то обхода. Требуется посчитать чёрную площадь внутри многоугольника.
Сначала я решил для случая, когда фигура — треугольник, а потом подумал, почему бы не триангулировать исходный многоугольник и не вызвать решение для треугольника $$$\left(n-2\right)$$$ раза?
Итак, всё упёрлось в триангуляцию. К сожалению, в статье выше многие моменты опускаются, есть всякие ошибки (к примеру, автор говорит, что вершина типа merge соединяется только с вершиной типа merge, хотя в примере, который он приводит, это не так, а в псевдокоде заифано на merge...). Плюс псевдокод, если его реализовывать так, как там написано, не проходит стресс-тестирование.
Я всё исправил и предоставляю свою реализацию. Итак, чтобы посчитать триангуляцию, надо вызвать следующую функцию:
std::vector<Triangle> triangulate(vpii points) {
makeCounterClockwise(points);
return triangulateFast(points);
}
Она вернёт вам вектор из треугольников, с которыми вы можете делать всё, что хотите. Будьте внимательны с тем, что я вырезаю вершины, которые лежат на одной прямой с соседними вершинами, ибо они бесполезны (т.е., если $$$P_{i-1}$$$, $$$P_{i}$$$ и $$$P_{i+1}$$$ лежат на одной прямой, то $$$P_i$$$ вырезается.
Исходный код (полное решение задачи H. Территориальное размежевание)
Исходный код
#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
#define isz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
void remin(auto &x, const auto &y) { if (x > y) x = y; }
void remax(auto &x, const auto &y) { if (x < y) x = y; }
using ll = long long;
using ld = long double;
using pii = std::pair<int,int>;
using vi = std::vector<int>;
using vpii = std::vector<pii>;
using vvpii = std::vector<vpii>;
using tiii = std::tuple<int,int,int>;
using vtiii = std::vector<tiii>;
#ifndef __GEOMA_HPP__
#define __GEOMA_HPP__
namespace algos {
namespace geoma {
template<typename A, typename B>
inline std::pair<A,B> operator-(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first - b.first, a.second - b.second);
}
template<typename A, typename B>
inline std::pair<A,B> operator+(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first + b.first, a.second + b.second);
}
template<typename T, typename A, typename B>
inline T cross(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return T(a.first) * b.second - T(a.second) * b.first;
}
template<typename T, typename A, typename B>
inline T doubleSquare(const std::pair<A,B> &a, const std::pair<A,B> &b, const std::pair<A,B> &c) {
T res = cross<T>(a-b, c-b);
if (res < 0)
res = -res;
return res;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A,B> &a) {
return T(a.first)*a.first + T(a.second)*a.second;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A, B> &a, const std::pair<A, B> &b) {
return dist2<T>(a-b);
}
using Triangle = std::tuple<pii,pii,pii>;
inline bool isPointInsideTriangle(Triangle tr, pii p) {
auto [A, B, C] = tr;
auto S = doubleSquare<ll>(A,B,C);
auto S1 = doubleSquare<ll>(p,B,C);
auto S2 = doubleSquare<ll>(A,p,C);
auto S3 = doubleSquare<ll>(A,B,p);
return S == S1 + S2 + S3;
}
inline ll orientedSquare(const std::vector<pii> &points) {
ll res = cross<ll>(points.back(), points.front());
for (int i = 0; i + 1 < isz(points); i++)
res += cross<ll>(points[i], points[i+1]);
return res;
}
inline vpii readPolygon() {
int n;
vpii points;
std::cin >> n;
points.resize(n);
for (auto &[x,y] : points)
std::cin >> x >> y;
return points;
}
inline void moveToPositive(vpii &points) {
int minX = INT_MAX, minY = INT_MAX;
for (auto &[x,y] : points) {
remin(minX, x);
remin(minY, y);
}
if (minX % 2 != 0) minX--;
if (minY % 2 != 0) minY--;
for (auto &[x,y] : points)
x -= minX, y -= minY;
}
inline void transpose(vpii &points) {
for (auto &[x,y] : points)
std::swap(x, y);
}
inline void makeCounterClockwise(vpii &points) {
ll sq = orientedSquare(points);
if (sq < 0) std::reverse(all(points));
}
inline int sign(auto x) { return x < 0 ? -1 : (x > 0 ? +1 : 0); }
const int LEFT = -1, RIGHT = +1;
inline int rotation(pii first, pii mid, pii last) {
// > 0 --> right
// < 0 --> left
return sign(cross<ll>(first-mid,last-mid));
}
inline bool isPointOnSegment(pii P, pii A, pii B) {
auto [x, y] = P;
auto [xmin,xmax] = std::minmax({A.first, B.first});
auto [ymin,ymax] = std::minmax({A.second, B.second});
if (x < xmin || x > xmax || y < ymin || y > ymax)
return false;
return cross<ll>(A-P,B-P) == 0;
}
const int NONE = 0, POINT = 1, INSIDE = 2;
inline int intersects(pii A, pii B, pii C, pii D) {
if (A == C || A == D || B == C || B == D)
return POINT;
if (isPointOnSegment(A,C,D) || isPointOnSegment(B,C,D) ||
isPointOnSegment(C,A,B) || isPointOnSegment(D,A,B))
return POINT;
auto [xmin1,xmax1] = std::minmax({A.first, B.first});
auto [ymin1,ymax1] = std::minmax({A.second, B.second});
auto [xmin2,xmax2] = std::minmax({C.first, D.first});
auto [ymin2,ymax2] = std::minmax({C.second, D.second});
if (xmax1 < xmin2 || xmin1 > xmax2)
return NONE;
if (ymax1 < ymin2 || ymin1 > ymax2)
return NONE;
auto rot1 = rotation(A,B,C);
auto rot2 = rotation(A,B,D);
if (sign(rot1) * sign(rot2) < 0)
return INSIDE;
return NONE;
}
inline bool isAlmostEqual(ld x, ld y) {
const ld EPS = 1e-10;
return std::abs(x-y) / (1 + std::abs(x)) < EPS;
}
inline bool angleLessThanPi(pii L, pii M, pii R) {
return rotation(L,M,R) == LEFT;
}
inline bool angleGreaterThanPi(pii L, pii M, pii R) {
return rotation(L,M,R) == RIGHT;
}
inline auto getPrev(const auto &vec, int i) {
i--;
if (i < 0) i = isz(vec)-1;
return vec[i];
}
inline auto getNext(const auto &vec, int i) {
i++;
if (i == isz(vec))
i = 0;
return vec[i];
}
inline bool isConvex(const vpii &P) {
for (int i = 0; i < isz(P); i++)
if (rotation(getPrev(P, i), P[i], getNext(P, i)) == RIGHT)
return false;
return true;
}
inline void removeSameLine(vpii &P) {
static std::vector<bool> keep;
keep.assign(isz(P), false);
for (int i = 0; i < isz(P); i++)
keep[i] = (rotation(getPrev(P,i), P[i], getNext(P,i)) != 0);
int p = 0;
for (int i = 0; i < isz(P); i++)
if (keep[i]) P[p++] = P[i];
P.resize(p);
}
} // namespace geoma
} // namespace algos
#endif // __GEOMA_HPP__
#ifndef __TRIANGULATE_HPP__
#define __TRIANGULATE_HPP__
namespace algos {
namespace triangulate {
using namespace geoma;
inline bool liesBelow(pii P1, pii P2) {
auto [x1,y1] = P1;
auto [x2,y2] = P2;
return (y1 < y2 || (y1 == y2 && x1 > x2));
}
inline bool liesAbove(pii P1, pii P2) {
auto [x1,y1] = P1;
auto [x2,y2] = P2;
return (y1 > y2 || (y1 == y2 && x1 < x2));
}
const int START = 0, SPLIT = 1, END = 2, MERGE = 3, REGULAR = 4;
inline auto getTypeOfVertex(const auto &P, int i) {
pii L = getPrev(P, i), M = P[i], R = getNext(P, i);
if (liesBelow(L,M) && liesBelow(R,M) && angleLessThanPi(L,M,R))
return START;
if (liesBelow(L,M) && liesBelow(R,M) && angleGreaterThanPi(L,M,R))
return SPLIT;
if (liesAbove(L,M) && liesAbove(R,M) && angleLessThanPi(L,M,R))
return END;
if (liesAbove(L,M) && liesAbove(R,M) && angleGreaterThanPi(L,M,R))
return MERGE;
return REGULAR;
}
inline std::string type2str(int t) {
switch(t) {
case START: return "start";
case SPLIT: return "split";
case END: return "end";
case MERGE: return "merge";
case REGULAR: return "regular";
default: return "none";
};
}
inline bool isMonotone(const vpii &P) {
if (isz(P) <= 3)
return true;
// we should start from max vertex
int imax = -1, ymax = INT_MIN;
for (int i = 0; i < isz(P); i++)
if (P[i].second >= ymax) {
imax = i;
ymax = P[i].second;
}
int len = 0;
int i = imax;
while (len <= isz(P) && P[i].second >= getNext(P, i).second) {
len++;
i++;
if (i == isz(P)) i = 0;
}
while (len <= isz(P) && P[i].second <= getNext(P, i).second) {
len++;
i++;
if (i == isz(P)) i = 0;
}
return len >= isz(P);
}
inline vpii findDiagonals(const vpii &P) {
static vtiii pq;
pq.clear();
for (int i = 0; i < isz(P); i++) {
auto [x,y] = P[i];
pq.emplace_back(y,-x,i);
}
std::sort(all(pq), std::greater<>());
int currY, currX;
auto setCurrY = [&](int newCurrY) { currY = newCurrY; };
auto setCurrX = [&](int newCurrX) { currX = newCurrX; };
static vi helper;
helper.assign(isz(P),-1);
auto getHelper = [&](int i) {
if (i < 0) i += isz(P);
assert(0 <= i && i < isz(P));
assert(helper[i] != -1);
return helper[i];
};
auto setHelper = [&](int i, int j) {
if (i < 0) i += isz(P);
assert(0 <= i && i < isz(P));
if (helper[i] == i)
helper[i] = j;
else if (helper[i] != -1) {
if (int k = helper[i]; P[k].second > P[j].second)
helper[i] = j;
} else
helper[i] = j;
};
auto getEdge = [&](int i) {
if (i < 0) i += isz(P);
if (i >= isz(P)) i -= isz(P);
auto L = P[i], R = getNext(P, i);
auto [ymin, ymax] = std::minmax({L.second, R.second});
return tiii(ymin, ymax, i);
};
auto getEdgeId = [&](int i) {
auto [y1,y2,j] = getEdge(i);
if (y1 == y2) return -1;
return j;
};
int cachedId = -1;
ld cachedX = -1;
auto isCached = [&](int id) { return id == cachedId; };
auto setCached = [&](int id, ld x) { cachedId = id; cachedX = x; };
auto clearCached = [&](){cachedId = -1; cachedX = -1;};
auto coordOfIntersection = [&](int i, int y) -> ld {
if (i < 0) i += isz(P);
if (isCached(i))
return cachedX;
assert(0 <= i && i < isz(P));
auto [x1,y1] = P[i];
auto [x2,y2] = getNext(P, i);
//ld B = ld(x2-x1)/(y2-y1);
//ld A = x1 - y1 * B;
// x1 - y1 * (x2-x1)/(y2-y1)
// x1 - y1 * (x2-x1)/(y2-y1) + y * (x2-x1)/(y2-y1)
// (x1 * (y2 - y1) + (y - y1) * (x2 - x1)) / (y2 - y1)
// (x1 * y2 + y * x2 - y * x1 - y1 * x2) / (y2 - y1)
assert(y2 != y1);
return ld((ld)x1 * y2 + (ld)y * x2 - (ld)y * x1 - (ld)y1 * x2) / (y2 - y1);
};
auto LessByIntersection = [&](const int &j1, const int &j2) {
ld x1 = coordOfIntersection(j1, currY);
ld x2 = coordOfIntersection(j2, currY);
if (x1 > x2) return false;
if (x1 < x2) return true;
return j1 < j2;
};
// y1 + t * (y2-y1) = y => t = (y - y1)/(y2 - y1)
// x1 + t * (x2-x1) = x => x1 + (y - y1) * (x2-x1) / (y2 - y1)
// coeff = (x2 - x1) / (y2 - y1)
// x1 + (y - y1) * coeff <= x
// x1 - y1 * coeff + y * coeff <= x
// A + y * B <= x -- need to sort by these values
std::set<int, decltype(LessByIntersection)> T(LessByIntersection);
auto insertEdgeInT = [&](int i) {
int ei = getEdgeId(i);
if (ei == -1) return;
auto y = P[i].second;
setCurrY(y);
setCached(ei, coordOfIntersection(ei, currY));
T.insert(ei);
clearCached();
};
auto deleteEdgeFromT = [&](int i) {
int j = getEdgeId(i);
if (j == -1) return;
auto it = T.find(j);
assert(it != T.end());
T.erase(it);
};
auto searchEdgeInT = [&](int i) {
assert(0 <= i && i < isz(P));
const auto [x, y] = P[i];
setCurrY(y); setCurrX(x);
int cand;
setCached(isz(P), x);
auto it = T.upper_bound(isz(P));
assert(it != T.begin());
cand = *std::prev(it);
clearCached();
return cand;
};
vpii diagonals;
auto addEdgeInTriangles = [&](int u, int v) {
if (getNext(P, u) != P[v] && getPrev(P, u) != P[v])
diagonals.emplace_back(u,v);
};
for (auto [y,x,i] : pq) {
x *= -1;
setCurrY(y);
switch (getTypeOfVertex(P,i)) {
case START: {
insertEdgeInT(i);
setHelper(i, i);
break;
}
case END: {
addEdgeInTriangles(i, getHelper(i-1));
deleteEdgeFromT(i-1);
break;
}
case SPLIT: {
int j = searchEdgeInT(i);
addEdgeInTriangles(i, getHelper(j));
setHelper(j, i);
insertEdgeInT(i);
setHelper(i, i);
break;
}
case MERGE: {
addEdgeInTriangles(i, getHelper(i-1));
deleteEdgeFromT(i-1);
int j = searchEdgeInT(i);
addEdgeInTriangles(i, getHelper(j));
setHelper(j, i);
break;
}
case REGULAR: {
if (liesAbove(getPrev(P,i), P[i]) && liesBelow(getNext(P,i), P[i])) {
addEdgeInTriangles(i, getHelper(i-1));
deleteEdgeFromT(i-1);
insertEdgeInT(i);
setHelper(i, i);
} else {
int j = searchEdgeInT(i);
addEdgeInTriangles(i, getHelper(j));
setHelper(j, i);
}
break;
}
};
}
return diagonals;
}
struct Node {
const vpii * parent = nullptr;
Node * next = nullptr;
Node * prev = nullptr;
int where{}, id{};
std::list<int> diagonals{};
void setNext(Node *n) {
next = n;
n->prev = this;
}
void setPrev(Node *n) {
prev = n;
n->next = this;
}
vpii to_vector() const {
vpii result(1, (*parent)[id]);
for (Node *it = next; it != this; it = it->next)
result.push_back((*parent)[it->id]);
return result;
}
pii to_point() const {
return (*parent)[id];
}
};
template<typename Func>
inline void splitByDiag(const vpii &P, const vpii &D, const Func &f) {
if (isz(P) <= 2)
return;
if (isz(P) == 3) {
assert(D.empty());
f(P);
return;
}
static std::vector<Node *> nodes;
nodes.clear();
for (int i = 0; i < isz(P); i++)
nodes.emplace_back(new Node{&P, nullptr, nullptr, 0, i, {}});
for (int i = 0; i+1 < isz(nodes); i++)
nodes[i]->setNext(nodes[i+1]);
nodes.front()->setPrev(nodes.back());
static std::vector<std::pair<Node*, Node*>> diagonals;
diagonals.clear();
for (const auto &[u,v] : D) {
assert(0 <= u && u < isz(nodes));
assert(0 <= v && v < isz(nodes));
int currDiag = isz(diagonals);
diagonals.push_back({nodes[u], nodes[v]});
nodes[u]->diagonals.push_back(currDiag);
nodes[v]->diagonals.push_back(currDiag);
}
int cntPoly = 1;
for (int id = 0; id < isz(diagonals); id++) {
// need to split by this diagonal
auto [startFi, startSe] = diagonals[id];
auto currFi = startFi, currSe = startSe;
while (currFi != startSe && currSe != startFi) {
currFi = currFi->next;
currSe = currSe->next;
}
// need to copy ends of diagonal:
nodes.push_back(new Node{&P, nullptr, nullptr, cntPoly, startFi->id, {}});
Node *copyFi = nodes.back();
nodes.push_back(new Node{&P, nullptr, nullptr, cntPoly, startSe->id, {}});
Node *copySe = nodes.back();
// choose a list with min len:
if (currFi != startSe) {
std::swap(copyFi, copySe);
std::swap(startFi, startSe);
std::swap(currFi, currSe);
}
assert(currFi == startSe);
// iterate over all vertices and updates diagonals
for (currFi = currFi->prev; currFi != startFi; currFi = currFi->prev) {
auto &dd = currFi->diagonals;
for (auto it = dd.begin(); it != dd.end(); ) {
// this diagonal still exists?
if (*it <= id) {
it = dd.erase(it);
continue;
}
auto &[lhs, rhs] = diagonals[*it];
if (lhs == startFi) (lhs = copyFi)->diagonals.push_back(*it);
if (rhs == startFi) (rhs = copyFi)->diagonals.push_back(*it);
if (lhs == startSe) (lhs = copySe)->diagonals.push_back(*it);
if (rhs == startSe) (rhs = copySe)->diagonals.push_back(*it);
it = std::next(it);
}
// set where this vertex lies
currFi->where = cntPoly;
}
// disconnect less half:
copyFi->setNext(currFi->next);
copySe->next = copyFi;
copySe->setNext(copyFi);
startSe->prev->setNext(copySe);
startFi->setNext(startSe);
cntPoly++;
}
static std::vector<Node *> anyNode;
anyNode.resize(cntPoly);
for (auto it : nodes)
anyNode[it->where] = it;
for (int i = 0; i < cntPoly; i++) {
vpii temp = anyNode[i]->to_vector();
removeSameLine(temp);
//assert(isMonotone(temp));
f(temp);
}
for (auto node : nodes)
delete node;
}
inline vvpii splitByDiag(const vpii &P, const vpii &D) {
vvpii result;
splitByDiag(P, D, [&](const vpii &PP){ result.push_back(PP); });
return result;
}
inline vpii TriangulateMonotonePolygon(const vpii &P) {
if (isz(P) <= 3)
return {};
// we need to create chains
int ymax = INT_MIN, ymin = INT_MAX, imax = -1;
int leftChainBegin = -1, leftChainEnd = -1;
for (int i = 0; i < isz(P); i++) {
if (P[i].second >= ymax) {
imax = i;
ymax = P[i].second;
}
if (P[i].second == ymax && getNext(P, i).second != ymax) {
leftChainBegin = i;
}
remin(ymin, P[i].second);
}
const int LEFT_CHAIN = 0, RIGHT_CHAIN = 1;
assert(imax != -1);
static vi chainType; chainType.assign(isz(P), -1);
static vi chainPos; chainPos.assign(isz(P), -1);
assert(leftChainBegin != -1);
leftChainEnd = leftChainBegin;
int pos = 0;
while (!(P[leftChainEnd].second == ymin && getPrev(P, leftChainEnd).second != ymin)) {
chainType[leftChainEnd] = LEFT_CHAIN;
chainPos[leftChainEnd] = pos++;
leftChainEnd++;
if (leftChainEnd == isz(P))
leftChainEnd = 0;
}
chainType[leftChainEnd] = LEFT_CHAIN;
if (chainPos[leftChainEnd] == -1)
chainPos[leftChainEnd] = pos++;
pos = 0;
for (int i = (leftChainBegin-1+isz(P)) % isz(P); chainType[i] == -1; ) {
chainType[i] = RIGHT_CHAIN;
chainPos[i] = pos++;
i--;
if (i < 0) i = isz(P)-1;
}
// {y-coord, chainId, posInChain}
auto cmp = [&](pii a, pii b) {
auto [ya, ida] = a;
auto [yb, idb] = b;
if (ya < yb) return true;
if (ya > yb) return false;
if (chainType[ida] != chainType[idb]) {
return chainType[ida] < chainType[idb];
}
return chainPos[ida] > chainPos[idb];
};
static vpii pq;
pq.clear();
for (int i = 0; i < isz(P); i++) {
auto [x,y] = P[i];
pq.push_back({y,i});
}
std::sort(pq.rbegin(), pq.rend(), cmp);
static vi S;
S.clear();
S.push_back(std::get<1>(pq[0]));
S.push_back(std::get<1>(pq[1]));
vpii diagonals;
auto addDiagonal = [&](int i, int j) {
if (getNext(P, i) != P[j] && getPrev(P, i) != P[j]) {
diagonals.emplace_back(i, j);
return true;
}
return false;
};
auto isSameChain = [&](int i, int j) {
return chainType[i] == chainType[j];
};
for (const auto &[y, id] : pq) {
if (!isSameChain(id,S.back())) {
int last = S.back();
while (isz(S) > 1) {
addDiagonal(id, S.back());
S.pop_back();
}
S = {last, id};
} else {
int want = LEFT, skip = RIGHT;
if (chainType[id] == LEFT_CHAIN)
std::swap(want, skip);
pii first = P[id];
while (isz(S) >= 2) {
pii mid = P[S[isz(S)-1]];
pii last = P[S[isz(S)-2]];
ll rot = rotation(first, mid, last);
if (rot == skip)
break;
if (rot != want) {
S.pop_back();
continue;
}
assert(rot == want);
addDiagonal(id, S[isz(S)-2]);
S.pop_back();
}
S.push_back(id);
}
}
return diagonals;
}
template<typename Func>
inline void triangulateFast(vpii inP, const Func &f) {
ll fullSquare = orientedSquare(inP);
removeSameLine(inP);
vpii D = findDiagonals(inP);
splitByDiag(inP, D, [&](vpii PP){
removeSameLine(PP);
vpii DD = TriangulateMonotonePolygon(PP);
auto ff = [&](const vpii &t){
fullSquare -= orientedSquare(t);
f(t);
};
splitByDiag(PP, DD, ff);
});
assert(fullSquare == 0);
}
inline std::vector<Triangle> triangulateFast(const vpii &inP) {
std::vector<Triangle> answ;
triangulateFast(inP, [&](const vpii &t){
assert(isz(t) <= 3);
if (isz(t) == 3) answ.emplace_back(t[0], t[1], t[2]);
});
return answ;
}
inline std::vector<Triangle> triangulateSlow(const vpii &points)
{
// slow triangulation in time O(n^2)
std::vector<Triangle> triangles;
std::list<pii> list(all(points));
while (isz(list) > 3) {
// search a vertex which we can cut off:
for (auto it = list.begin(); it != list.end(); ) {
auto next = (std::next(it) == list.end() ? list.begin() : std::next(it));
auto prev = (it == list.begin() ? std::prev(list.end()) : std::prev(it));
ll rot = rotation(*prev, *it, *next);
if (rot == 0) {
it = list.erase(it);
} else if (rot < 0) {
// check this vertex:
bool ok = true;
for (auto p : list) {
bool bad = isPointInsideTriangle({*prev, *it, *next}, p);
bad = bad && !isPointOnSegment(p, *prev, *it);
bad = bad && !isPointOnSegment(p, *it, *next);
bad = bad && !isPointOnSegment(p, *prev, *next);
if (bad) {
ok = false;
break;
}
}
if (ok) {
triangles.emplace_back(*prev, *it, *next);
it = list.erase(it);
} else {
it++;
}
} else {
it++;
}
}
}
if (isz(list) == 3)
triangles.emplace_back(*list.begin(), *(++list.begin()), *(++++list.begin()));
else assert(isz(list) < 3);
return triangles;
}
inline std::vector<Triangle> triangulateConvex(const vpii &points) {
// points is convex polygon
// take any vertex and make triangles
std::vector<Triangle> triangles;
for (int i = 1; i + 1 < isz(points); i++)
triangles.emplace_back(points[0], points[i], points[i+1]);
return triangles;
}
inline std::vector<Triangle> triangulate(vpii points) {
makeCounterClockwise(points);
return triangulateFast(points);
}
} // namespace triangulate
} // namespace algos
#endif // __TRIANGULATE_HPP__
using namespace algos::triangulate;
using namespace algos::geoma;
pii pickPointOutside(pii A, pii B) {
pii C1(A.first, B.second);
pii C2(B.first, A.second);
pii res = cross<ll>(A-C1, B-C1) < 0 ? C1 : C2;
return res;
}
std::pair<ld,ld> rectangle(int xmin, int ymin, int xmax, int ymax) {
ll dx = xmax - xmin, dy = ymax - ymin;
ld fi = (dx + 1) / 2 * dy;
ld se = dx / 2 * dy;
if (xmin % 2 != 0) std::swap(fi,se);
return {fi, se};
}
ll sum(ll n) {
// n^2 - (n-1)^2 + (n-2)^2 - (n-3)^2 + (n-4)^2 + ...
return n * (n+1) / 2;
}
std::pair<ld, ld> solveBaseCaseLD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
ld c = y / (2.0L * x);
std::pair<ld, ld> res{c * sum(x), c * sum(x-1)};
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseRU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseRD(xmax-xmin,ymax-ymin);
if (xmin % 2 != 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> triangle(pii A, pii B, pii C) {
// B - corner
// AC - diagonal
if (doubleSquare<ll>(A,B,C) == 0)
return {0, 0};
if (A.first > C.first)
std::swap(A, C);
ll xmin = std::min({A.first, B.first, C.first});
ll xmax = std::max({A.first, B.first, C.first});
ll ymin = std::min({A.second, B.second, C.second});
ll ymax = std::max({A.second, B.second, C.second});
pii LD(xmin,ymin), LU(xmin,ymax), RD(xmax,ymin), RU(xmax,ymax);
std::pair<ld,ld> res;
if (A.second > C.second) {
if (B == LD) res = solveBaseCaseLD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseRU(xmin,ymin,xmax,ymax);
} else {
if (B == RD) res = solveBaseCaseRD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseLU(xmin,ymin,xmax,ymax);
}
return res;
}
std::pair<ld,ld> twoRectangles(pii A, pii B, pii C) {
if (A.first > B.first) std::swap(A, B);
if (A.first > C.first) std::swap(A, C);
if (B.first > C.first) std::swap(B, C);
assert(A.first <= B.first && B.first <= C.first);
auto [ymin, ymax] = std::minmax({A.second,C.second});
if (B.second >= ymax || B.second <= ymin || B.first == A.first || B.first == C.first) {
remin(ymin, B.second);
remax(ymax, B.second);
return rectangle(A.first, ymin, C.first, ymax);
}
std::pair<ld,ld> rect = rectangle(A.first, ymin, C.first, ymax);
// above or below than diagonal?
if (rotation(A,C,B) <= 0) { // above
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(A.first, B.second, B.first, ymax);
else
rect = rect - rectangle(B.first, B.second, C.first, ymax);
} else { // below
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(B.first, ymin, C.first, B.second);
else
rect = rect - rectangle(A.first, ymin, B.first, B.second);
}
return rect;
}
std::pair<ld,ld> calcTriangle(pii A, pii B, pii C) {
if (cross<ll>(A-B,C-B) > 0)
std::swap(A, C);
auto temp = twoRectangles(A,B,C);
auto P1 = pickPointOutside(A,B);
auto P2 = pickPointOutside(B,C);
auto P3 = pickPointOutside(C,A);
if (!(A.first == B.first || A.second == B.second))
temp = temp - triangle(A,P1,B);
if (!(B.first == C.first || B.second == C.second))
temp = temp - triangle(B,P2,C);
if (!(C.first == A.first || C.second == A.second))
temp = temp - triangle(C,P3,A);
return temp;
}
ld solveFast(vpii points) {
moveToPositive(points);
transpose(points);
makeCounterClockwise(points);
ld answ = 0;
triangulateFast(points, [&](const vpii &tr){
assert(isz(tr) <= 3);
if (isz(tr) == 3)
answ += calcTriangle(tr[0],tr[1],tr[2]).first;
});
return answ;
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(0);
std::cout << std::fixed << std::setprecision(12);
auto points = readPolygon();
std::cout << solveFast(points) << std::endl;
}
Это решение работает от $$$280$$$ до $$$296$$$ мс в системе Яндекс.Контест при $$$n=10^5$$$.
UPD. Обновил исходный код, удалив все неиспользуемые функции и закомментированные строки.
UPD 2. В комментариях на русском языке orz предложил решение оригинальной задачи за $$$O(N)$$$. Давайте применим тот же трюк, который делается при вычислении обычной площади многоугольника. Зафиксируем точку $$$P=(0,0)$$$ и посчитаем ориентированные площади треугольников $$$P A_i A_{i+1}$$$ для каждой стороны: если при обходе точек в данном порядке был левый поворот, то прибавим к ответу чёрную площадь внутри этого треугольника, иначе вычтем её из ответа.
Linear solution of problem H
#pragma GCC optimize("Ofast")
#include <bits/stdc++.h>
#define isz(x) (int)(x).size()
#define all(x) (x).begin(), (x).end()
void remin(auto &x, const auto &y) { if (x > y) x = y; }
void remax(auto &x, const auto &y) { if (x < y) x = y; }
using ll = long long;
using ld = long double;
using pii = std::pair<int,int>;
using vi = std::vector<int>;
using vpii = std::vector<pii>;
using vvpii = std::vector<vpii>;
using tiii = std::tuple<int,int,int>;
using vtiii = std::vector<tiii>;
#ifndef __GEOMA_HPP__
#define __GEOMA_HPP__
namespace algos {
namespace geoma {
template<typename A, typename B>
inline std::pair<A,B> operator-(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first - b.first, a.second - b.second);
}
template<typename A, typename B>
inline std::pair<A,B> operator+(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return std::pair<A,B>(a.first + b.first, a.second + b.second);
}
template<typename T, typename A, typename B>
inline T cross(const std::pair<A,B> &a, const std::pair<A,B> &b) {
return T(a.first) * b.second - T(a.second) * b.first;
}
template<typename T, typename A, typename B>
inline T doubleSquare(const std::pair<A,B> &a, const std::pair<A,B> &b, const std::pair<A,B> &c) {
T res = cross<T>(a-b, c-b);
if (res < 0)
res = -res;
return res;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A,B> &a) {
return T(a.first)*a.first + T(a.second)*a.second;
}
template<typename T, typename A, typename B>
inline T dist2(const std::pair<A, B> &a, const std::pair<A, B> &b) {
return dist2<T>(a-b);
}
inline ll orientedSquare(const std::vector<pii> &points) {
ll res = cross<ll>(points.back(), points.front());
for (int i = 0; i + 1 < isz(points); i++)
res += cross<ll>(points[i], points[i+1]);
return res;
}
inline vpii readPolygon() {
int n;
vpii points;
std::cin >> n;
points.resize(n);
for (auto &[x,y] : points)
std::cin >> x >> y;
return points;
}
inline void moveToPositive(vpii &points) {
int minX = INT_MAX, minY = INT_MAX;
for (auto &[x,y] : points) {
remin(minX, x);
remin(minY, y);
}
if (minX % 2 != 0) minX--;
if (minY % 2 != 0) minY--;
for (auto &[x,y] : points)
x -= minX, y -= minY;
}
inline void transpose(vpii &points) {
for (auto &[x,y] : points)
std::swap(x, y);
}
inline void makeCounterClockwise(vpii &points) {
ll sq = orientedSquare(points);
if (sq < 0) std::reverse(all(points));
}
inline int sign(auto x) { return x < 0 ? -1 : (x > 0 ? +1 : 0); }
const int LEFT = -1, RIGHT = +1;
inline int rotation(pii first, pii mid, pii last) {
// > 0 --> right
// < 0 --> left
return sign(cross<ll>(first-mid,last-mid));
}
} // geoma
} // algos
#endif // __GEOMA_HPP__
using namespace algos::geoma;
pii pickPointOutside(pii A, pii B) {
pii C1(A.first, B.second);
pii C2(B.first, A.second);
pii res = cross<ll>(A-C1, B-C1) < 0 ? C1 : C2;
return res;
}
std::pair<ld,ld> rectangle(int xmin, int ymin, int xmax, int ymax) {
ll dx = xmax - xmin, dy = ymax - ymin;
ld fi = (dx + 1) / 2 * dy;
ld se = dx / 2 * dy;
if (xmin % 2 != 0) std::swap(fi,se);
return {fi, se};
}
ll sum(ll n) {
// n^2 - (n-1)^2 + (n-2)^2 - (n-3)^2 + (n-4)^2 + ...
return n * (n+1) / 2;
}
std::pair<ld, ld> solveBaseCaseLD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
ld c = y / (2.0L * x);
std::pair<ld, ld> res{c * sum(x), c * sum(x-1)};
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseRU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll x, ll y) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseLD(x, y);
if (x % 2 == 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll x, ll y) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLD(x, y);
return res;
}
std::pair<ld, ld> solveBaseCaseLU(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// ****
// ***
// **
// *
auto res = solveBaseCaseLU(xmax-xmin, ymax-ymin);
if (xmin % 2 != 0) std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> solveBaseCaseRD(ll xmin, ll ymin, ll xmax, ll ymax) {
// triangle
// *
// **
// ***
// ****
auto res = solveBaseCaseRD(xmax-xmin,ymax-ymin);
if (xmin % 2 != 0)
std::swap(res.first, res.second);
return res;
}
std::pair<ld, ld> triangle(pii A, pii B, pii C) {
// B - corner
// AC - diagonal
if (doubleSquare<ll>(A,B,C) == 0)
return {0, 0};
if (A.first > C.first)
std::swap(A, C);
ll xmin = std::min({A.first, B.first, C.first});
ll xmax = std::max({A.first, B.first, C.first});
ll ymin = std::min({A.second, B.second, C.second});
ll ymax = std::max({A.second, B.second, C.second});
pii LD(xmin,ymin), LU(xmin,ymax), RD(xmax,ymin), RU(xmax,ymax);
std::pair<ld,ld> res;
if (A.second > C.second) {
if (B == LD) res = solveBaseCaseLD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseRU(xmin,ymin,xmax,ymax);
} else {
if (B == RD) res = solveBaseCaseRD(xmin,ymin,xmax,ymax);
else res = solveBaseCaseLU(xmin,ymin,xmax,ymax);
}
return res;
}
std::pair<ld,ld> twoRectangles(pii A, pii B, pii C) {
if (A.first > B.first) std::swap(A, B);
if (A.first > C.first) std::swap(A, C);
if (B.first > C.first) std::swap(B, C);
assert(A.first <= B.first && B.first <= C.first);
auto [ymin, ymax] = std::minmax({A.second,C.second});
if (B.second >= ymax || B.second <= ymin || B.first == A.first || B.first == C.first) {
remin(ymin, B.second);
remax(ymax, B.second);
return rectangle(A.first, ymin, C.first, ymax);
}
std::pair<ld,ld> rect = rectangle(A.first, ymin, C.first, ymax);
// above or below than diagonal?
if (rotation(A,C,B) <= 0) { // above
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(A.first, B.second, B.first, ymax);
else
rect = rect - rectangle(B.first, B.second, C.first, ymax);
} else { // below
// is A below than C?
if (A.second <= C.second)
rect = rect - rectangle(B.first, ymin, C.first, B.second);
else
rect = rect - rectangle(A.first, ymin, B.first, B.second);
}
return rect;
}
std::pair<ld,ld> calcTriangle(pii A, pii B, pii C) {
if (cross<ll>(A-B,C-B) > 0)
std::swap(A, C);
auto temp = twoRectangles(A,B,C);
auto P1 = pickPointOutside(A,B);
auto P2 = pickPointOutside(B,C);
auto P3 = pickPointOutside(C,A);
if (!(A.first == B.first || A.second == B.second))
temp = temp - triangle(A,P1,B);
if (!(B.first == C.first || B.second == C.second))
temp = temp - triangle(B,P2,C);
if (!(C.first == A.first || C.second == A.second))
temp = temp - triangle(C,P3,A);
return temp;
}
ld solveFast(vpii points) {
moveToPositive(points);
transpose(points);
makeCounterClockwise(points);
ld answ = 0;
pii A(0,0);
for (int i = 0; i < isz(points); i++) {
pii B = points[i], C = points[(i+1)%isz(points)];
ld curr = calcTriangle(A, B, C).first;
if (rotation(A,B,C) == LEFT) answ += curr;
else answ -= curr;
}
return answ;
}
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(0);
std::cout << std::fixed << std::setprecision(12);
auto points = readPolygon();
std::cout << solveFast(points) << std::endl;
}
UPD 3. k1r1t0 предложил использовать вместо треугольников трапеции, основания которых параллельны осям Ox или Oy. Я пока не понял, как складывать и вычитать трапеции, поэтому исходника нет.
