Hello, codeforces!
Today, based on the article Triangulation of polygon from a book "Computational Geometry: Algorithms and applications", I implemented an algorithm of triangulation of any polygon without self intersections (non-convex polygons too) in C++.
$$$\text{ }$$$ $$$\text{ }$$$
Why it was needed for me? Probably, you know about a competition RuCode 7.0. In this year, qualification round contained next problem H: you are given an endless plane, painted in black and white stripes with height 1 and infinity width. These stripes alternate: black, white, black, white... You are given a polygon on this plane as list of $$$n\text{ }\left(3\le n \le 100000\right)$$$ 2D points with integer coordinates $$$x_i, y_i \left(-10^9 \le x_i, y_i \le 10^9\right)$$$ in some traversal order. Need to calculate a total square inside this polygon which is painted in black color.
Firstly, I solved for any triangle ($$$n=3$$$). Secondly, I thought about the triangulation of polygon and accumulating answers for each triangle.
So, the heaviest part of solution is the triangulation of polygon. I implemented this, stress-tested and fixed all of incomprehensible and dubious places of pseudo code that can be found in above book. For example, the author of this book states that vertex which has a type "merge" should be connected with another "merge" vertex, and this is used in pseudo-code, but there is an example, where it is not true. Stress-testing shows all of corner cases, which are not covered by the pseudo code.
I fixed all of mistakes which the stress-testing has found. You can just call this function:
std::vector<Triangle> triangulate(vpii points) {
makeCounterClockwise(points);
return triangulateFast(points);
}
The output of this function is a vector of triangles. Be careful, because I'm erasing non-meaningful vertices like if $$$P_{i-1}$$$, $$$P_{i}$$$ and $$$P_{i+1}$$$ lie on the same line than erase $$$P_i$$$.
Source code of solution for given problem
#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;
}
This solution works in $$$296$$$ ms in Yandex.Contest when $$$n=10^5$$$.
UPD. updated the source code by removing all unused things and commented lines
UPD 2. orz supposed $$$O(N)$$$ solution of original problem H in russian comments for this blog. Lets use same trick which is used in an algorithm of calculation of square of polygon. We can fix a point $$$P=(0,0)$$$ and calculate the answer for each triangle $$$P A_i A_{i+1}$$$ (for each side of polygon): if the order of traversak of current triangle was clockwise, then we can take this answer with a sign "minus", otherwise we can take this answer with a sign "plus".
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 supposed to use trapezoids instead of triangles (in russian comments for this blog). Seems that for trapezoids a square which is painted in black can be calculated simplier.
