J: Cover the Polygon with Your Disk
問題
円と凸多角形が与えられる。凸多角形の頂点は10個以下。それぞれ自由に平行移動・回転移動させて、共通部分の面積の最大値を出力せよ。
入力は整数で座標は0~100、円の半径は1~100。
解法
円の中心が凸多角形の内部にある場合だけ考えればよい。幾何の証明は詳しくないが、凸多角形なので図示すれば明らか。
円の中心を凸多角形の内部(x, y)で移動させ、その共通部分の面積を2引数関数f(x, y)とすると、fは定義域に対してUnimodal （もしかしたら凸）になりそうなので三分探索したい。2次元の三分探索はx方向・y方向の2つの探索を入れ子にする。
共通部分の面積パートについて。円の中心が多角形の内部にあるとすると、中心と多角形の頂点に線を伸ばしてピザカットした扇型と三角形の共通部分を求めれば良い。
その後、円と線分の交点を求めて共通部分を更に小さな三角形と扇型に分解すればよい。
コード
三分探索は均等に3分割する必要はなく、中心に近い二点を取れば収束が早くなるはず。この実装では黄金比分割して計算回数を減らそうとしている。
#include<complex> #include<stdio.h> #include<iostream> #include<vector> #include<algorithm> #include<string> #include<string.h> #include<cmath> using namespace std; typedef long long LL; typedef vector<int> VI; #define REP(i,n) for(int i=0, i##_len=(n); i<i##_len; ++i) #define EACH(i,c) for(__typeof((c).begin()) i=(c).begin(),i##_end=(c).end();i!=i##_end;++i) #define eprintf(s...) fprintf(stderr, s) template<class T> inline void amin(T &x, const T &y) { if (y<x) x=y; } template<class T> inline void amax(T &x, const T &y) { if (x<y) x=y; } const double EPS = 1e-8; // 0 ~~ [-EPS, EPS] const double INF = 1e12; const double PI = acos((double)-1); typedef complex<double> Point; const Point I(0, 1); // abs(P a) ::= sqrt(a.x^2 + a.y^2); // norm(P a) ::= a.x^2 + a.y^2; // Point polar(rho, theta=0) int sign(double x) { if (x < -EPS) return -1; if (EPS < x) return 1; return 0; } int cmp(double x, double y) { return sign(x - y); } namespace std { bool operator<(const Point &x, const Point &y) { return real(x) != real(y)? real(x) < real(y): imag(x) < imag(y); } } double cross(const Point &x, const Point &y) { return imag(conj(x)*y); } double dot(const Point &x, const Point &y) { return real(conj(x)*y); } Point unit_vec(const Point &a) { return a / abs(a); } struct Line { Point p[2]; Line() {} Line(const Line &l) { p[0] = l[0]; p[1] = l[1]; } Line(const Point &a, const Point &b) { p[0] = a; p[1] = b; } Point& operator[](int i) { return p[i]; } const Point& operator[](int i) const { return p[i]; } }; Point unit_vec(const Line &l) { return unit_vec(l[1] - l[0]); } // normal line Point normal(const Line &l) { return unit_vec((l[1] - l[0]) * I); } struct Circle { Point p; double r; Circle(){} Circle(const Point &p, double (r)) : p(p), r(r) {} }; string toStr(double x) { static char buf[60]; sprintf(buf, "%.6f", x); return buf; } string toStr(const Point &p) { return string("(") + toStr(p.real()) + ", " + toStr(p.imag()) + ")"; } string toStr(const Line &l) { return string("(") + toStr(l[0]) + " -- " + toStr(l[1]) + ")"; } string toStr(const Circle &c) { return string("(Circle ") + toStr(c.p) + ", " + toStr(c.r) + ")"; } Point projection(const Line &l, const Point &p) { double t = dot(p-l[0], l[0]-l[1]) / norm(l[0]-l[1]); return l[0]+t*(l[0]-l[1]); } vector<Point> crossPoint(const Line &l, const Circle &c) { vector<Point> ret; Point p = projection(l, c.p); double d = abs(p - c.p); if (c.r < d) return ret; double s = sqrt(c.r * c.r - d * d); Point v = unit_vec(l) * s; ret.push_back(p - v); ret.push_back(p + v); return ret; } struct ByArg { bool operator()(const Point &x, const Point &y) const { return arg(x) < arg(y); } } byArg; int N; double R; Point Z[21]; double area(Point a, Point b) { if (abs(a) <= R && abs(b) <= R) return cross(a, b) / 2; double arg_ = arg(a); a *= polar(1.0, -arg_); b *= polar(1.0, -arg_); vector<Point> c = crossPoint(Line(a, b), Circle(Point(), R)); sort(c.begin(), c.end(), byArg); c.push_back(b); double ret = 0; Point p = a; EACH (e, c) { if (cmp(abs(a - *e) + abs(b - *e), abs(a - b)) == 0) { if (cmp(abs(p), R) <= 0 && cmp(abs(*e), R) <= 0) { ret += cross(p, *e) / 2; } else { ret += R * R * (arg(*e) - arg(p)) / 2; } p = *e; } } return ret; } double area() { Point p = Z[N-1]; double ret = 0; REP (i, N) { ret += area(p, Z[i]); p = Z[i]; } return ret; } const double GR = (sqrt(5) - 1.0) / 2.0; // GR == 0.618 double argmax(double f(double), double lo, double hi) { double x = hi - GR * (hi - lo), y = lo + GR * (hi - lo); double vx = f(x), vy = f(y); for (int i=0; i<60; i++) { if (vx < vy) { lo = x; x = y; vx = vy; y = lo + GR * (hi - lo); vy = f(y); } else { hi = y; y = x; vy = vx; x = hi - GR * (hi - lo); vx = f(x); } } return (x+y)/2.0; } double y_shift(double dy) { REP (i, N) Z[i].imag() += dy; double ret = area(); REP (i, N) Z[i].imag() -= dy; return ret; } double x_shift(double dx) { REP (i, N) Z[i].real() += dx; Point p = Z[N-1]; double bound_D = INF, bound_U = -INF; REP (i, N) { if (cmp(p.real(), 0) <= 0 && cmp(0, Z[i].real()) <= 0) { double y = (Z[i].imag() - p.imag()) * (-p.real()) / (Z[i].real() - p.real()) + p.imag(); if (cmp(Z[i].real(), p.real()) == 0) y = Z[i].imag(); amax(bound_U, -y); } if (cmp(0, p.real()) <= 0 && cmp(Z[i].real(), 0) <= 0) { double y = (Z[i].imag() - p.imag()) * (-p.real()) / (Z[i].real() - p.real()) + p.imag(); if (cmp(Z[i].real(), p.real()) == 0) y = Z[i].imag(); amin(bound_D, -y); } p = Z[i]; } double dy = argmax(y_shift, bound_D, bound_U); REP (i, N) Z[i].imag() += dy; double ret = area(); REP (i, N) Z[i].real() -= dx; return ret; } int main() { scanf("%d", &N); { int r; scanf("%d", &r); R = r; } REP (i, N) { int x, y; scanf("%d%d", &x, &y); Z[i] = Point(x, y); } double bound_L = INF, bound_R = -INF; REP (i, N) { amin(bound_L, -Z[i].real()); amax(bound_R, -Z[i].real()); } double dx = argmax(x_shift, bound_L, bound_R); REP (i, N) Z[i].real() += dx; printf("%.9f\n", area()); return 0; }
感想
円の中心が多角形の外側にある場合や、多角形が凸ではない場合、多角形を符号付き面積を持つ三角形に分解することで共通部分の面積を求めることができるらしい。
ただ、円の中心が外側にあると三分探索での定義域を求めるのがやや難しい。ちなみに凸とは限らない問題では三分探索ができない。
幾何は得意ではないが昨年と比べ簡単で良かった。
国内予選の幾何のほうが10倍難しく、昨年の幾何のほうが100倍難しい（私見）。AOJを見ると昨年の幾何は3000~5500-byte 程度で解かれているので、なにか良い方法があったのかもしれない。
