FOSSIL (상, p. 486)
본 글은 과거에 작성한 풀이 기록이다.#
이후 동일 문제를 보다 간결하고 명확한 방식으로 다시 정리하였다.
본 글의 구현 방식은 혼자서 생각해낸 것이 많아 다소 장황하며, 불필요한 로직이 구현되어 있고, 구현 방식이 깔끔하지 않다.
따라서 최신 글을 참고하는 것이 좋으며, 본 글은 지우기 아까워 남겨놓는다.
최신 글 보러 가기 ☞ FOSSIL (상, p. 486) - 다시 풀어봄
문제#
어떻게 풀었나?#
두 Convex Hull이 겹치는 구간도 Convex Hull이다. 이 새로운 Convex Hull에서 남북 방향의 최대 길이를 구하면 되는데, 남북 방향의 최대 길이는 왼쪽에서부터 오른쪽으로 갈 수록 점점 커졌다가 줄어들 수 밖에 없다. 그래서 그 길이를 삼분 탐색(Ternary Search로 구하면 된다.
상세 방법은 아래와 같다.
1. Convex Hull을 12시 방향의 점부터 반시계 방향 순서대로 정렬한다.#
정렬하는 이유는 Convex Hull을 점이 아닌 선분기준으로 봐야하기 때문이다.
두 Convex Hull 간의 교차점을 구하기 위해선, 점이 아닌 선분을 기준으로 반복문을 돌아야 한다. 그래서 Convex Hull을 이루는 꼭지점들을 반시계 방향 순서로 정렬해야 [i]번째 꼭지점과 [i+1]번째 꼭지점이 선분을 이룰 수 있다.
또 어떤 점이 한 Convex Hull 안에 포함되는 경우를 구할 때에도 Convex Hull의 꼭지점들은 정렬되어 있어야 한다. 왜냐하면, (다른 방법이 있는지는 모르겠지만) Convex Hull 안에 포함되는지 여부를 확인하기위한 점을 기준으로 Convex Hull들의 꼭지점들을 잇는 선분을 그리고 해당 선분들을 반시계 방향 순으로 돌면서 각도가 $\pi$를 초과하는지 확인해야 하기 때문이다. 각도가 $\pi$를 초과하는 경우 그 점은 Convex Hull안에 포함될 수 없다.
이 외에도 마지막에 구한 새로운 Convex Hull을 삼분 탐색할 때에도 반시계 방향 정렬이 필요하다.
참고로 위에서 언급한 반시계 방향은 모두 시계 방향으로도 가능하다. 글 작성 편의상 한쪽 방향만 언급했다.
2. 두 Convex Hull의 교차점을 모두 구한다.#
교차점을 구하는 방식이 중요한데, 벡터의 외적을 활용한다.
1) 먼저 교차점의 존재 여부부터 파악해야 한다.#
2차원에서 두 벡터 $\vec{u}=(x_1,y_1),\; \vec{v}=(x_2,y_2)$의 외적은 다음과 같은 스칼라 값으로 정의한다.
$$ \vec{u} \times \vec{v} = x_1y_2 - y_1x_2 $$그리고 이 결과값의 부호가 방향을 의미한다.
> 0: 반시계(CCW)< 0: 시계(CW)= 0: 같은 직선(colinear)
참고로 이 결과값의 크기(절대값)는 두 벡터가 이루는 평행사변형의 넓이를 의미한다.
아래 그림을 보면, $\overline{AB}$와 $\overline{CD}$는 교차한다.
그러면 $\overrightarrow{AB} \times \overrightarrow{AC} \;$와 $\overrightarrow{AB} \times \overrightarrow{AD} \;$는 부호가 다를 수밖에 없다.
2) 교차점이 존재 여부를 파악했으면, 교차점($I$) 을 구한다.#
$\overline{AB}$위의 모든 점은 아래와 같이 벡터와 파라미터를 사용해서 나타낼 수 있다.
$$ P(t) = A + t(B-A),\quad 0\le t \le 1 $$$\overline{CD}$도 마찬가지로 아래와 같이 나타낼 수 있다.
$$ Q(u) = C + u(D-C), \quad 0 \le u \le 1 $$그러므로 교차점 $I$는 아래의 식을 풀면 된다.
$$ A + t(B-A) = C+u(D-C) $$$$ \Rightarrow t(B-A) = (C-A) + u(D-C) $$양변에 $(D-C)$로 외적을 취한다.
$$ \Rightarrow t(B-A)\times (D-C) = (C-A)\times(D-C) $$왜냐면, $u(D-C) \times (D-C) = 0$ 이기 때문이다. (자기 자신과의 외적은 항상 0)
$$ \therefore t = \frac{(C-A)\times(D-C)}{(B-A)\times (D-C)} $$여기서 우리는 분모가 0이 된다는 것은 $\overline{AB} \parallel \overline{CD}$인데, 위의 방식으로 교차점을 구하기 전에 교차점이 있는지 이미 확인한 상태이므로 분모가 0이 되는 상황은 발생하지 않는다.
이제 위의 t값을 P(t)에 대입하면 교차점 $I$를 구할 수 있다.
3. 한 Convex Hull이 다른 Convex Hull의 꼭지점을 포함하는 경우도 모두 구한다.#
1번에서 Convex Hull의 꼭지점을 반시계 방향으로 정렬하는 이유에 대해서 설명할 때에도 얘기했지만, 좀 더 그림으로 자세히 설명하면 아래와 같다.
아래 그림의 점P가 Convex Hull ABCDE 내부의 점이려면, 점P에서 A부터 E까지 반시계 방향 순서로 선분들을 이어줄 때, 이웃한 선분들의 사이각 $\theta$값이 $\pi$를 초과해서는 안된다.
이 점을 이용하여 Convex Hull 안에 들어가는 다른 Convex Hull의 꼭지점들을 모두 구해준다.
4. 2, 3번에서 구한 점들이 겹치는 구간을 의미하는 새로운 Convex Hull이 된다.#
만약 해당 점이 2개 이하라면 Convex Hull을 구성할 수 없으므로, 남북 방향의 최대 길이는 0이 된다.
5. 새로운 Convex Hull을 12시 방향의 점부터 반시계 방향 순서대로 정렬한다.#
6. 새로운 Convex Hull에서 남북 방향의 최대 길이를 삼분 탐색으로 구한다.#
정답 코드#
#include <bits/stdc++.h>
using namespace std;
using pdd = pair<double, double>;
inline double PI() {
static double pi = acos(-1);
return pi;
}
vector<pdd> v, w;
int GetIdx12clockCoord(vector<pdd>& g) {
double mx = -101;
int mxi = -1;
for (int i = 0; i < g.size(); i++) {
if (mx < g[i].second) {
mx = g[i].second;
mxi = i;
}
}
return mxi;
}
inline double GetTheta(double x, double y) {
return fmod(atan2(y, x) + 2 * PI(), 2 * PI());
}
void SortByConvexHullOrder(vector<pdd>& g) {
if (g.empty()) return;
//12시 방향 제일 윗 좌표를 고른다.
int idx = GetIdx12clockCoord(g);
if (idx != 0) swap(g[0], g[idx]);
double gijun = GetTheta(0, 1);
for (int i = 0; i < g.size() - 1; i++) {
double mn_nt = 2 * PI();
double mnj = -1;
for (int j = i + 1; j < g.size(); j++) {
double nx = g[j].first - g[i].first;
double ny = g[j].second - g[i].second;
double theta = GetTheta(nx, ny);
double nt = fmod(theta - gijun + 4 * PI(), 2 * PI());
if (mn_nt > nt) {
mn_nt = nt;
mnj = j;
}
}
if (i + 1 != mnj) swap(g[i + 1], g[mnj]);
gijun += mn_nt;
}
}
void GetIncCo(vector<pdd>& in, vector<pdd>& out, vector<pdd>& ret) {
//in.size()-1까지 반복문을 도는 것은 in[0]가 in vector의 마지막에 한번 더 추가되었기 때문이다.
//이 처리는 out vector에도 동일하게 적용되어있다.
for (int i = 0; i < in.size() - 1; i++) {
double prv_theta = 0;
bool isInc = true;
for (int j = 0; j < out.size(); j++) {
double nx = out[j].first - in[i].first;
double ny = out[j].second - in[i].second;
double theta = GetTheta(nx, ny);
if (j != 0) {
double nt = fmod(theta - prv_theta + 2 * PI(), 2 * PI());
if (nt > PI()) {
isInc = false;
break;
}
}
prv_theta = theta;
}
if (isInc) ret.push_back(in[i]);
}
}
inline pdd VectorMinus(pdd x, pdd y) {
return pdd(x.first - y.first, x.second - y.second);
}
inline pdd VectorPlus(pdd x, pdd y) {
return pdd(x.first + y.first, x.second + y.second);
}
inline pdd VectorMul(double t, pdd x) {
return pdd(t * x.first, t * x.second);
}
inline int Sign(double x) {
return (x > -1e-8) - (x < 1e-8);
}
inline double Cross(pdd x, pdd y) {
return x.first * y.second - x.second * y.first;
}
inline int Ccw(pdd x, pdd y) {
return Sign(Cross(x, y));
}
void GetIntersection(int vi, int vj, int wi, int wj, vector<pdd>& cvh) {
pdd a = VectorMinus(v[vj], v[vi]);
pdd b = VectorMinus(v[vj], w[wi]);
pdd c = VectorMinus(v[vj], w[wj]);
pdd d = VectorMinus(w[wj], w[wi]);
pdd e = VectorMinus(w[wj], v[vi]);
pdd f = VectorMinus(w[wj], v[vj]);
//교점이 있는 경우
if (Ccw(a, b) * Ccw(a, c) <= 0 && Ccw(d, e) * Ccw(d, f) <= 0) {
double r = Cross(VectorMinus(v[vj], v[vi]), VectorMinus(w[wj], w[wi]));
if (fabs(r) < 1e-8) return; //교점이 없음.
double q = Cross(VectorMinus(w[wi], v[vi]), VectorMinus(w[wj], w[wi]));
double t = q / r;
pdd insec = VectorPlus(v[vi], VectorMul(t, VectorMinus(v[vj], v[vi])));
cvh.push_back(insec);
}
}
double F(vector<pdd>& g, double x) {
static vector<double> ans;
ans.clear();
for (int i = 0; i < g.size() - 1; i++) {
if ((g[i].first <= x && x <= g[i + 1].first) ||
(g[i + 1].first <= x && x <= g[i].first)) {
if (fabs(g[i].first - g[i + 1].first) < 1e-8) { //==0
ans.push_back(g[i].second);
ans.push_back(g[i + 1].second);
}
else if (fabs(g[i].second - g[i + 1].second) < 1e-8) {
ans.push_back(g[i].second);
}
else {
double a = (g[i].second - g[i + 1].second) / (g[i].first - g[i + 1].first);
ans.push_back(a * (x - g[i].first) + g[i].second);
}
}
}
if (ans.empty()) return 0.0;
sort(ans.begin(), ans.end());
return ans.back() - ans[0];
}
int main() {
ios::sync_with_stdio(false); cin.tie(0);
int tc; cin >> tc;
while (tc--) {
int n, m;
cin >> n >> m;
v.clear(); w.clear();
for (int i = 0; i < n; i++) {
double x, y;
cin >> x >> y;
v.emplace_back(x, y);
}
for (int i = 0; i < m; i++) {
double x, y;
cin >> x >> y;
w.emplace_back(x, y);
}
//처음 고른 좌표부터 Convex Hull을 그리는 좌표 순서로 정렬한다.
SortByConvexHullOrder(v);
SortByConvexHullOrder(w);
//GetIntersection()에서 Convex Hull의 선분을 한바퀴 돌아야 하기 때문에
//가장 첫번째 점인 v[0],w[0]를 맨 마지막에도 한번 더 넣어준다.
v.push_back(v[0]);
w.push_back(w[0]);
static vector<pdd> cvh;
cvh.clear();
//두 ConvexHull의 모든 교점 구하기
for (int i = 0; i < v.size() - 1; i++) {
for (int j = 0; j < w.size() - 1; j++) {
GetIntersection(i, i + 1, j, j + 1, cvh);
}
}
//한쪽 ConvexHull의 점이 다른 한쪽의 ConvexHull 안으로 들어가는 경우 cvh에 추가.
GetIncCo(v, w, cvh);
GetIncCo(w, v, cvh);
SortByConvexHullOrder(cvh);
double ans = 0.0;
if (cvh.size() > 2) { //ternary search
double lo = 101, hi = -1.0;
for (int i = 0; i < cvh.size(); i++) {
lo = min(lo, cvh[i].first);
hi = max(hi, cvh[i].first);
}
cvh.push_back(cvh[0]);
for (int z = 0; z < 100; z++) {
double lmid = (lo * 2 + hi) / 3;
double rmid = (lo + hi * 2) / 3;
double lres = F(cvh, lmid);
double rres = F(cvh, rmid);
if (lres < rres) lo = lmid;
else hi = rmid;
}
ans = F(cvh, lo);
}
printf("%.10lf\n", ans);
}
return 0;
}주의할 점#
아래와 같은 경우도 고려되었는지 확인해야 한다.
좌측의 그림은 두 Convex Hull간의 겹치는 부분에 두 Convex Hull의 꼭지점이 하나도 포함되어 있지 않은 경우도 있다는 것을 말한다.
우측의 그림은 한 Convex Hull이 다른 Convex Hull을 완전히 포함하는 경우를 말한다.
추가 Test Case와 Output#
//Input
2
3 3
0 0 6 2 0 4
2 2 4 0 4 4
4 4
20 20 80 20 80 80 20 80
40 40 60 40 60 60 40 60
//Output
2.00000000
20.00000000Comment#
처음으로 풀어본 제대로된 기하문제가 아니었나 싶다. 전에 풀어봤던 MINASTRITH 문제는 원을 직선으로 바꿔 풀어서 벡터를 다룰 일이 없었으니 진짜 기하문제를 다뤄봤다고 얘기하긴 힘들어 보였는데, 이 문제는 정말 순수한(?) 기하 문제라 말할 수 있을 것 같다.
문제를 보고 아이디어는 바로 떠올랐다. 두 컨벡스 헐의 겹치는 부분이 새로운 컨벡스 헐을 보여주고, 그 새로운 컨벡스 헐의 남북방향 최대 길이는 삼분 탐색으로 찾으면 되는게 보였다. 문제는 새로운 컨벡스 헐을 코드로 구해내는 것이었다. 이 부분은 일부 챗지의 도움을 받았다. 위에서 두 컨벡스 헐의 교차점을 벡터로 구하는 방식을 챗지가 알려줬다. 그 외 나머지 부분은 모두 책을 안보고 풀어서 뿌듯하다. 컨벡스 헐 안에 포함되는 점을 판별하는 것도 고민하다가 떠올린 방법인데, 책은 어떻게 설명하고 있는지 모르겠다.
이제 책 코드를 보러 가야겠다. 비슷하게 풀었으려나..?
책 풀이#
책을 보고 왔다. 완전히 잘못 풀었다. 아, 아니다. 잘못 푼건 아니지만, 내가 너무 어렵게 풀었다. 이 문제는 기하 문제 푸는 방법들을 사용하지 않고도 풀 수 있는 문제였다. 생각해보니 내가 보고 있던 책의 챕터도 기하쪽이 아니라 수치 해석쪽이었다.
위에서 나는 겹치는 구간의 도형을 완벽히 구한 다음, 거기에서 삼분탐색을 통해 가장 긴 남북방향의 길이를 찾았다. 그런데 책에서는 두 ConvexHull의 윗껍질과 아랫껍질로 나눠서 x값을 대입하면 윗껍질에선 가장 작은 y값을 찾고, 아랫껍질에선 가장 큰 y값을 찾아서 두 차이를 x좌표에 해당하는 겹치는 구간의 길이로 봤다.
메인 아이디어는 이렇고, 그럼 예외처리가 복잡하지 않을까 싶은데.. 생각해보면 위에서 구한 구간의 길이가 음수가 되는 경우 겹치는 경우가 아니며, 최대 내각이 $\pi$를 넘지 않는 Convex Hull이기 때문에 x좌표계에 수직인 선분의 경우 해당 선분의 x좌표 지점에서 남북방향 최대 길이가 답이 되는 경우는 절대 발생할 수 없다. 따라서 이런 $x=a$ 꼴의 선분의 경우 윗껍질 or 아랫껍질중 어디에 넣을지 고민하지 않아도 된다.
삼분탐색을 할 때에도 lo, hi값 세팅할 때 굳이 겹치는 지점을 명확히 구해서 시작할 필요가 없었다. 두 Convex Hull의 가장 작은 x값중 max값을 선택해 lo를 지정하면 됐고, hi는 두 Convex Hull의 가장 큰 x값중 min값을 선택해 hi를 지정하면 됐다. 이렇게 해서 lmid나 rmid쪽이 겹치는 구간이 아니라면, lres나 rres쪽이 음수가 되면서 음수가 된 쪽의 lo or hi가 정답쪽으로 좁혀지게 된다. 그래서 이와 관련된 예외처리를 구구절절 코드로 만들어줄 필요가 없다.
책 풀이 방식으로 풀어본 정답 코드#
#include <bits/stdc++.h>
using namespace std;
using pdd = pair<double, double>;
using ppddpdd = pair<pdd, pdd>;
vector<pdd> v, w;
vector<ppddpdd> upper, lower;
bool IsBetween(double x1, double x2, double x) {
return (x1 <= x && x <= x2) || (x2 <= x && x <= x1);
}
double At(ppddpdd& line, double x) {
double a = (line.second.second - line.first.second) / (line.second.first - line.first.first);
return a * (x - line.first.first) + line.first.second;
}
double GetUpperY(vector<ppddpdd>& upper, double x) {
double mn_y = 100;
for (int i = 0; i < upper.size(); i++) {
if (IsBetween(upper[i].first.first, upper[i].second.first, x)) {
mn_y = min(mn_y, At(upper[i], x));
}
}
return mn_y;
}
double GetLowerY(vector<ppddpdd>& lower, double x) {
double mx_y = 0;
for (int i = 0; i < lower.size(); i++) {
if (IsBetween(lower[i].first.first, lower[i].second.first, x)) {
mx_y = max(mx_y, At(lower[i], x));
}
}
return mx_y;
}
double GetMinX(vector<pdd>& g) {
double mn = 100;
for (int i = 0; i < g.size() - 1; i++) {
mn = min(mn, g[i].first);
}
return mn;
}
double GetMaxX(vector<pdd>& g) {
double mx = 0;
for (int i = 0; i < g.size() - 1; i++) {
mx = max(mx, g[i].first);
}
return mx;
}
int main() {
ios::sync_with_stdio(false); cin.tie(0);
int tc; cin >> tc;
while (tc--) {
int n, m;
cin >> n >> m;
v.clear(); w.clear();
lower.clear(); upper.clear();
double x, y;
for (int i = 0; i < n; i++) {
cin >> x >> y;
v.emplace_back(x, y);
}
for (int i = 0; i < m; i++) {
cin >> x >> y;
w.emplace_back(x, y);
}
v.push_back(v[0]);
w.push_back(w[0]);
for (int i = 1; i < v.size(); i++) {
if (v[i].first > v[i - 1].first) lower.emplace_back(v[i - 1], v[i]);
else if (v[i].first < v[i - 1].first) upper.emplace_back(v[i], v[i - 1]);
}
for (int i = 1; i < w.size(); i++) {
if (w[i].first > w[i - 1].first) lower.emplace_back(w[i - 1], w[i]);
else if (w[i].first < w[i - 1].first) upper.emplace_back(w[i], w[i - 1]);
}
sort(lower.begin(), lower.end());
sort(upper.begin(), upper.end());
double lo = max(GetMinX(v), GetMinX(w));
double hi = min(GetMaxX(v), GetMaxX(w));
for (int z = 0; z < 100; z++) {
double lmid = (lo * 2 + hi) / 3;
double rmid = (lo + hi * 2) / 3;
double lres = GetUpperY(upper, lmid) - GetLowerY(lower, lmid);
double rres = GetUpperY(upper, rmid) - GetLowerY(lower, rmid);
if (lres < rres) lo = lmid;
else hi = rmid;
}
double ans_x = lo;
double ans_len = GetUpperY(upper, ans_x) - GetLowerY(lower, ans_x);
if (ans_len < 0) ans_len = 0.0;
printf("%.8lf\n", ans_len);
}
return 0;
}