2019-2020 ICPC Asia Taipei-Hsinchu Regional Contest - L Largest Quadrilateral
Attachments - 2019-2020 ICPC Asia Taipei-Hsinchu Regional Contest - Codeforces
概要
N点与えられる。そのうちの適切な順番で4つの点を選び、作れる最大の四角形の面積を求めよ。
制約： 4 <= N <= 4096, 0 <= X_i, Y_i <= 10^9
ACするまで
1. なるべく大きい四角形を考えると、凸包に含まれてる点から選ばないと行けなさそう。
2. わかりません。
3. 凸包は、時には円のように考えてみるとうまくいきやすい！
4. 作る四角形の対角線のうちの1本を固定するとする。固定したうえでの面積最大化を考えてみる。
5. その状態で、対角線の一端を固定、もう一端を少しずらしてみると、面積が最大化されるような四角形の残りの2つの頂点の動き方は、広義単調であると気付く。->変則しゃくとり(命名、僕)でいけそう 例外処理忘れずに(すべて1本の直線上とか、凸包上の点が3つしかないとか)
6. AC!幾何ライブラリを持っておくと、2次元計算が楽になるのでいいですね。
考察
まず、点の集合のうち、自明に凸包上にある点が候補になるとわかる。いろいろ面倒なのは考えずに、ここでは、凸包上の点の数が4つ以上存在してるとして考えていく。
凸包の上からいかにして4点を選ぶか？
という問題に対して、凸包では次のように考えるとうまくいきやすい。
Point
凸包上での動き方は、円の上での動き方を似た感じに適用できる場合が多い。
今回の場合は、この考えでうまくいく。つまり、次のような問題を考えてみる。
円上に何点かあり、それらが作る四角形の最大面積を求めよ。
大学受験でのテクニックではあるが、円に内接する四角形は、対角線を1本固定すると、それを底辺とした三角形2つそれぞれの面積最大化を独立に考えればよい。今回もこの考えでいく。
atcoder.jp
紙の最後で言及してるしゃくとり法は、この問題の尺取法と同じ。次の記事で詳しく書かれている。
sen-comp.hatenablog.com
ACソースコード
#include<iostream> #include<algorithm> #include<cmath> #include<vector> using namespace std; typedef long long ll; namespace v2d { struct Vec2 { //2次元ベクトルのクラス double x, y; Vec2(double _x, double _y) { x = _x, y = _y; } Vec2() { x = 0, y = 0; } Vec2 operator+() const { return *this; } Vec2 operator-() const { return{ -x, -y }; } Vec2 operator+ (const Vec2& b) const { return{ x + b.x, y + b.y }; } Vec2 operator- (const Vec2& b) const { return{ x - b.x, y - b.y }; } Vec2 operator* (const double b) const { return{ x * b, y * b }; } Vec2 operator/ (const double b) const { return{ x / b, y / b }; } Vec2 operator+= (const Vec2& b) { x += b.x, y += b.y; return *this; } Vec2 operator-= (const Vec2& b) { x -= b.x, y -= b.y; return *this; } Vec2 operator*= (const double b) { x *= b, y *= b; return *this; } Vec2 operator/= (const double b) { x /= b, y /= b; return *this; } bool operator== (const Vec2& b) { return b.x == x && b.y == y; } double lengthSquare() { return (x * x + y * y); } double length() { return std::sqrt(lengthSquare()); } double dot(const Vec2& b) { return x * b.x + y * b.y; } double cross(const Vec2& b) { //Generally, cross product is vector, but in 2D, cross product is also scalar. return x * b.y - y * b.x; } double distanceFrom(const Vec2& b) { return std::sqrt((x - b.x) * (x - b.x) + (y - b.y) * (y - b.y)); } Vec2 normalized() { return{ x / length(), y / length() }; } bool isZero() { return x == 0.0 && y == 0.0; } }; inline Vec2 operator*(double a, const Vec2& b) { return{ b.x * a, b.y * a }; } template <class Char> inline std::basic_ostream<Char>& operator <<(std::basic_ostream<Char>& os, const Vec2& v) { return os << Char('(') << v.x << Char(',') << v.y << Char(')'); } template <class Char> inline std::basic_istream<Char>& operator >> (std::basic_istream<Char>& is, Vec2& v) { return is >> v.x >> v.y; } } using namespace v2d; namespace convex_h { bool comp_x(const v2d::Vec2& p, const v2d::Vec2& q) { if (p.x != q.x)return p.x < q.x; return p.y < q.y; } bool comp_y(const v2d::Vec2& p, const v2d::Vec2& q) { if (p.y != q.y)return p.y < q.y; return p.x < q.x; } struct convex_hull { const double EPS = 1e-10; vector<int> ans; vector<v2d::Vec2> candidate; vector<pair<v2d::Vec2, int>> taiou; void calc(bool is_cover = true) { for (int i = 0; i < candidate.size(); i++) { taiou.push_back(make_pair(candidate[i], i)); } sort(taiou.begin(), taiou.end(), [](pair<v2d::Vec2, int> &a, pair<v2d::Vec2, int> &b) { return comp_x(a.first, b.first); }); int k = 0; ans.resize(2 * taiou.size()); for (int i = 0; i < taiou.size(); i++) { while (k > 1 && (is_cover && ((candidate[ans[k - 1]] - candidate[ans[k - 2]]). cross(taiou[i].first - candidate[ans[k - 1]]) <= 0) || (!is_cover && (candidate[ans[k - 1]] - candidate[ans[k - 2]]). cross(taiou[i].first - candidate[ans[k - 1]]) < 0)))k--; ans[k++] = taiou[i].second; } for (int i = taiou.size() - 2, t = k; i >= 0; i--) { while (k > t && (is_cover && (candidate[ans[k - 1]] - candidate[ans[k - 2]]). cross(taiou[i].first - candidate[ans[k - 1]]) <= 0) || ((!is_cover && (candidate[ans[k - 1]] - candidate[ans[k - 2]]). cross(taiou[i].first - candidate[ans[k - 1]]) < 0)))k--; ans[k++] = taiou[i].second; } ans.resize(k - 1); } double farthest_dis() { //Using effective algorithm, //this function can calculate farthest distance of all point in convex hull int O(n). //This algorithm name is Rptating Calipers, according to 蟻本. double res = 0; int n = ans.size(); if (ans.size() == 2) { return candidate[ans[0]].distanceFrom(candidate[ans[1]]); } int i = 0, j = 0; for (int k = 0; k < n; k++) { if (!comp_x(candidate[ans[i]], candidate[ans[k]]))i = k; if (comp_x(candidate[ans[j]], candidate[ans[k]]))j = k; } int si = i, sj = j; while (i != sj || j != si) { res = max(res, candidate[ans[i]].distanceFrom(candidate[ans[j]])); int nxti = (i + 1) % n, nxtj = (j + 1) % n; if ((candidate[ans[nxti]] - candidate[ans[i]]).cross(candidate[ans[nxtj]] - candidate[ans[j]]) < 0) { i = (i + 1) % n; } else { j = (j + 1) % n; } } return res; } double hull_distance() { double ret = candidate[ans[0]].distanceFrom(candidate[ans.back()]); for (int i = 1; i < ans.size(); i++) { ret += candidate[ans[i]].distanceFrom(candidate[ans[i - 1]]); } return ret; } }; } ll S_2times(Vec2 a, Vec2 b) { return abs(a.x * b.y - a.y * b.x); } ll min3(ll a, ll b, ll c) { return min(a, min(b, c)); } int main() { int T; cin >> T; while (T--) { int N; cin >> N; vector<Vec2> P(N); for (int i = 0; i < N; i++) { cin >> P[i]; } convex_h::convex_hull ch; ch.candidate = P, ch.calc(false); //P[ch.ans[i]] if (ch.ans.size() <= 2) { cout << 0 << endl; continue; } if (ch.ans.size() == 3) { ll ans = 0; Vec2 a = P[ch.ans[0]], b = P[ch.ans[1]], c = P[ch.ans[2]]; for (int i = 0; i < N; i++) { if (a == P[i] || b == P[i] || c == P[i])continue; ans = max(ans, S_2times(b - a, c - a) - min3( S_2times(P[i] - a, c - a), S_2times(P[i] - a, b - a), S_2times(P[i] - b, P[i] - c) ) ); } if (ans % 2 == 0)cout << ans / 2 << endl; else cout << ans / 2 << ".5" << endl; continue; } ll ans = -1145141919810; for (int i = 0; i < ch.ans.size(); i++) { int s1 = i, s2 = i + 2; ll maxs1 = 0, maxs2 = 0; for (int j = i + 2; j < ch.ans.size() + i; j++) { while (s1 + 1 < j && maxs1 <= S_2times(P[ch.ans[j % ch.ans.size()]] - P[ch.ans[i]], P[ch.ans[j % ch.ans.size()]] - P[ch.ans[(s1 + 1) % ch.ans.size()]])) { maxs1 = S_2times(P[ch.ans[j % ch.ans.size()]] - P[ch.ans[i]], P[ch.ans[j % ch.ans.size()]] - P[ch.ans[(s1 + 1) % ch.ans.size()]]); s1++; } while (s2 + 1 < ch.ans.size() + i && maxs2 <= S_2times(P[ch.ans[j % ch.ans.size()]] - P[ch.ans[i]], P[ch.ans[j % ch.ans.size()]] - P[ch.ans[(s2 + 1) % ch.ans.size()]])) { maxs2 = S_2times(P[ch.ans[j % ch.ans.size()]] - P[ch.ans[i]], P[ch.ans[j % ch.ans.size()]] - P[ch.ans[(s2 + 1) % ch.ans.size()]]); s2++; } ans = max(ans, maxs1 + maxs2); if (s2 == j)s2++; else { maxs1 = S_2times(P[ch.ans[(j + 1) % ch.ans.size()]] - P[ch.ans[i]], P[ch.ans[(j + 1) % ch.ans.size()]] - P[ch.ans[s1 % ch.ans.size()]]); maxs2 = S_2times(P[ch.ans[(j + 1) % ch.ans.size()]] - P[ch.ans[i]], P[ch.ans[(j + 1) % ch.ans.size()]] - P[ch.ans[s2 % ch.ans.size()]]); } } } if(ans % 2 == 0)cout << ans / 2 << endl; else cout << ans / 2 << ".5" << endl; } return 0; }
GYMの問題なので、普通の人は人の提出コードが見れない。なので、ソースコードは長いけど貼った。
2次元座標クラスは、cross()が2次元版外積に相当し、dot()がない積に相当。
凸包クラスのcalc()の引数は、falseならば、凸包に含まれる点で必ず凸包が曲がるようになる。
一言
このRegionalセット、解法の文献が乏しすぎる。日本語資料はじゅぴろさんの実況っぽいやつしかないし、中国語資料はコード、ぺ！ｗしかない。(ついでに言うと中国語資料のコードは理解できなかった)
解法のアイディアをくれたHyadoくん、ありがとう！
まとめ
Point
凸包上での動き方は、円の上での動き方を似た感じに適用できる場合が多い。
