Jan 6, 2024
Rant on Incorrect ECNA 2023 C (Computational Geometry)
So I was studying this ECNA 2023 C problem: Convex Hull Extension. And it just so turns out that none of the official solutions are correct…
Why Official Solutions are not Correct
To answer this question, we must first understand the official solution to the problem. After some analysis of the problem, we can see that the problem is actually asking us how many integer points ther are in the extension area of each edge, formed by extending the adjacent edges until they intersect. For example, for this shaded hexagon below, the extension area of edge \( BD \) would be \( \triangle BDP \).
Now this is all good and correct, but then the official solution falls short on this claim:
// ...
/**
* Since the bound on x,y coordinates is |x|,|y|<=1000, each triangle can
* cover a range of coordinates of size up to ~4 million. Calculating the
* x coordinate range for each y coordinate is can be done in constant time.
* Where N is the number of points and K is the bound on x,y coordinates,
* this approach has complexity O(N*K^2). For N=50 and K=1000, this should
* be feasible.
**/
// ...
To understand why this is not correct, let us consider these two segments: \( (0, 0)-(1000, 999) \) and \( (-1, -1000)-(1000, 0) \). Observe that the first segment has slope \( \frac{999}{1000} \) and the second segment has slope \( \frac{1000}{1001} \), so the difference in slope is around \( 10^{-6} \). Because the difference in intercept is around \( 10^3 \), two segments won't intersect until we extend them to around \( 10^9 \) in both x and y coordinate. This polygon below is built based on this principle:
5
998 997
0 0
-250 -500
-1 -1000
1000 0
And all official solutions just TLE because the coordinate range extends to \( 10^9 \). Note that some official solutions try to optimize by using Pick's theorem when the intersection point is an integer point. By ensuring that this does not happen, we can make all official solutions TLE.
Correct Solution
The correct solution here is to use Pick's theorem in a smart way. Notably, we divide \( \triangle BDP \) into quadrilateral \( BDD'B' \) and \( \triangle B'D'P \) such that integer points \( B' \) and \( D' \) are as close to \( P \) as possible. The exact value can be found in the graph below.
Observe that the number of integer points in quadrilateral \( BDD'B' \) can be found with Pick's theorem. And the coordinate span in \( \triangle B'D'P \) is less than \( 2\Delta x \), which makes finding the number of integer points in it solvable by brute force.
My solution is available below:
#include <bits/stdc++.h>
#define double __float128
#define cd const double &
const double EPS = 1E-9, PI = acos(-1);
int sgn(cd x) { return x < -EPS ? -1 : x > EPS; }
int cmp(cd x, cd y) { return sgn(x - y); }
double sqr(cd x) { return x * x; }
double msqrt(cd x) { return sgn(x) <= 0 ? 0 : sqrt(x); }
double abs(double x) { return x < 0 ? -x : x; }
#define cp const point &
struct point {
double x, y;
explicit point(cd x = 0, cd y = 0) : x(x), y(y) {}
int dim() const { return sgn(y) == 0 ? sgn(x) > 0 : sgn(y) > 0; }
point unit() const {
double l = msqrt(x * x + y * y);
return point(x / l, y / l);
}
point rot90() const { return point(-y, x); }
point _rot90() const { return point(y, -x); }
point rot(cd t) const {
double c = cos(t), s = sin(t);
return point(x * c - y * s, x * s + y * c);
}
};
bool operator==(cp a, cp b) { return cmp(a.x, b.x) == 0 && cmp(a.y, b.y) == 0; }
bool operator!=(cp a, cp b) { return cmp(a.x, b.x) != 0 || cmp(a.y, b.y) != 0; }
bool operator<(cp a, cp b) {
return cmp(a.x, b.x) == 0 ? cmp(a.y, b.y) < 0 : cmp(a.x, b.x) < 0;
}
point operator-(cp a) { return point(-a.x, -a.y); }
point operator+(cp a, cp b) { return point(a.x + b.x, a.y + b.y); }
point operator-(cp a, cp b) { return point(a.x - b.x, a.y - b.y); }
point operator*(cp a, cd b) { return point(a.x * b, a.y * b); }
point operator/(cp a, cd b) { return point(a.x / b, a.y / b); }
double dot(cp a, cp b) { return a.x * b.x + a.y * b.y; }
double det(cp a, cp b) { return a.x * b.y - a.y * b.x; }
double dis2(cp a, cp b = point()) { return sqr(a.x - b.x) + sqr(a.y - b.y); }
double dis(cp a, cp b = point()) { return msqrt(dis2(a, b)); }
#define cl const line &
struct line {
point s, t;
explicit line(cp s = point(), cp t = point()) : s(s), t(t) {}
};
bool point_on_segment(cp a, cl b) {
return sgn(det(a - b.s, b.t - b.s)) == 0 && sgn(dot(b.s - a, b.t - a)) <= 0;
}
bool two_side(cp a, cp b, cl c) {
return sgn(det(a - c.s, c.t - c.s)) * sgn(det(b - c.s, c.t - c.s)) < 0;
}
bool intersect_judgment(cl a, cl b) {
if (point_on_segment(b.s, a) || point_on_segment(b.t, a))
return true;
if (point_on_segment(a.s, b) || point_on_segment(a.t, b))
return true;
return two_side(a.s, a.t, b) && two_side(b.s, b.t, a);
}
point line_intersect(cl a, cl b) {
double s1 = det(a.t - a.s, b.s - a.s), s2 = det(a.t - a.s, b.t - a.s);
return (b.s * s2 - b.t * s1) / (s2 - s1);
}
double point_to_line(cp a, cl b) {
return abs(det(b.t - b.s, a - b.s)) / dis(b.s, b.t);
}
point project_to_line(cp a, cl b) {
return b.s + (b.t - b.s) * (dot(a - b.s, b.t - b.s) / dis2(b.t, b.s));
}
double point_to_segment(cp a, cl b) {
if (sgn(dot(b.s - a, b.t - b.s) * dot(b.t - a, b.t - b.s)) <= 0)
return abs(det(b.t - b.s, a - b.s)) / dis(b.s, b.t);
return std::min(dis(a, b.s), dis(a, b.t));
}
bool in_polygon(cp p, const std::vector<point> &po) {
int n = (int)po.size(), counter = 0;
for (int i = 0; i < n; ++i) {
point a = po[i], b = po[(i + 1) % n];
// Modify the next line if necessary.
if (point_on_segment(p, line(a, b)))
return true;
int x = sgn(det(p - a, b - a)), y = sgn(a.y - p.y), z = sgn(b.y - p.y);
if (x > 0 && y <= 0 && z > 0)
counter++;
if (x < 0 && z <= 0 && y > 0)
counter--;
}
return counter != 0;
}
double polygon_area(const std::vector<point> &a) {
double ans = 0.0;
for (int i = 0; i < (int)a.size(); ++i)
ans += det(a[i], a[(i + 1) % a.size()]) / 2.0;
return ans;
}
int N;
std::vector<point> hull;
long long point_on_boundary(point p) {
long long x = abs(p.x), y = abs(p.y);
return std::gcd(x, y) + 1;
}
long long point_on_boundary_poly(const std::vector<point> &p) {
long long ans = 0;
for (int i = 0; i < (int)p.size(); ++i) {
ans += point_on_boundary(p[(i + 1) % (int)p.size()] - p[i]) - 1;
}
return ans;
}
int main() {
scanf("%d", &N);
for (int i = 0; i < N; ++i) {
long long x, y;
scanf("%lld%lld", &x, &y);
hull.push_back(point(x, y));
}
bool is_inf = false;
long long ans = 0;
for (int i = 0; i < N; ++i) {
line l1 = line(hull[i], hull[(i + 1) % N]);
line l2 = line(hull[(i + 3) % N], hull[(i + 2) % N]);
if (sgn(det(l1.t - l1.s, l2.t - l2.s) >= 0)) {
std::vector<point> poly = {l1.t, l1.t + (l1.t - l1.s) * 2500 * 2500,
l2.t + (l2.t - l2.s) * 2500 * 2500,
l2.t};
long long point_in_poly =
polygon_area(poly) + 1 - point_on_boundary_poly(poly) / 2;
is_inf |= point_in_poly > 0;
continue;
}
point p = line_intersect(l1, l2);
double k1 = floor(dis(p - l1.t) / dis(l1.t - l1.s));
double k2 = floor(dis(p - l2.t) / dis(l2.t - l2.s));
std::vector<point> poly = {l1.t, l1.t + (l1.t - l1.s) * k1,
l2.t + (l2.t - l2.s) * k2, l2.t};
long long point_in_poly =
polygon_area(poly) + 1 - point_on_boundary_poly(poly) / 2;
if (p != poly[1] && p != poly[2])
point_in_poly +=
std::max(0ll, point_on_boundary(poly[2] - poly[1]) - 2);
std::vector<point> tri = {l1.t + (l1.t - l1.s) * k1,
l2.t + (l2.t - l2.s) * k2, p};
double min_x = std::min(std::min(tri[0].x, tri[1].x), tri[2].x);
double max_x = std::max(std::max(tri[0].x, tri[1].x), tri[2].x);
long long point_in_tri = 0;
for (long long i = ceil(min_x) + (!cmp(min_x, ceil(min_x)));
i <= floor(max_x) - (!cmp(max_x, floor(max_x))); ++i) {
std::vector<double> inter;
for (int j = 0; j < 3; ++j) {
point p1 = tri[j], p2 = tri[(j + 1) % 3];
if (cmp(std::min(p1.x, p2.x), i) <= 0 &&
cmp(i, std::max(p1.x, p2.x)) <= 0) {
inter.push_back(
line_intersect(line(p1, p2),
line(point(i, 0), point(i, 1)))
.y);
}
}
std::sort(inter.begin(), inter.end());
point_in_tri += std::max(
(long long)((floor(inter.back()) -
(!cmp(inter.back(), floor(inter.back())))) -
(ceil(inter.front()) +
(!cmp(inter.front(), ceil(inter.front())))) +
1),
0ll);
}
ans += point_in_poly + point_in_tri;
}
if (is_inf)
puts("infinitely many");
else
printf("%lld\n", ans);
}
Rant
I think ultimately this issue did not cause too much of an impact since this problem was very difficult to implement by itself. Almost no team got it in contest. Still, it is kind of frustrating to see that a much easier solution (i.e. the official one) could have gotten AC but the idea got rejected in the contest for thinking something that the judge didn't think of.
Also, our best player here at Purdue (@eggag32) spent two hours on this problem, and we did not qualify for NAC this year, so…
