#include using namespace std; using ll = long long; using ld = long double; // -------------------------------------------------------- #define FOR(i,l,r) for (int i = (l); i < (r); ++i) #define RFOR(i,l,r) for (int i = (r)-1; (l) <= i; --i) #define REP(i,n) FOR(i,0,n) #define RREP(i,n) RFOR(i,0,n) #define ALL(c) (c).begin(), (c).end() #define RALL(c) (c).rbegin(), (c).rend() #define SORT(c) sort(ALL(c)) #define RSORT(c) sort(RALL(c)) #define MIN(c) *min_element(ALL(c)) #define MAX(c) *max_element(ALL(c)) #define SUM(c) accumulate(ALL(c), 0LL) #define COUNT(c,v) count(ALL(c),(v)) #define SZ(c) ((ll)(c).size()) #define BIT(b,i) (((b)>>(i)) & 1) #define PCNT(b) __builtin_popcountll(b) #define P0(i) (((i) & 1) == 0) #define P1(i) (((i) & 1) == 1) #ifdef _LOCAL #define debug_bar cerr << "--------------------\n"; #define debug(x) cerr << "l." << __LINE__ << " : " << #x << " = " << (x) << '\n' #define debug_pair(x) cerr << "l." << __LINE__ << " : " << #x << " = (" << x.first << "," << x.second << ")\n"; template void debug_line(const vector& ans, int l, int r, int L = 0) { cerr << "l." << L << " :"; for (int i = l; i < r; i++) { cerr << ' ' << ans[i]; } cerr << '\n'; } #else #define cerr if (false) cerr #define debug_bar #define debug(x) #define debug_pair(x) template void debug_line([[maybe_unused]] const vector& ans, [[maybe_unused]] int l, [[maybe_unused]] int r, [[maybe_unused]] int L = 0) {} #endif template void input(T&... a) { (cin >> ... >> a); } void print() { cout << '\n'; } template void print(const T& a) { cout << a << '\n'; } template void print(const T& a, const Ts&... b) { cout << a; (cout << ... << (cout << ' ', b)); cout << '\n'; } template void cout_line(const vector& ans, int l, int r) { for (int i = l; i < r; i++) { if (i != l) { cout << ' '; } cout << ans[i]; } cout << '\n'; } template bool chmin(T& a, const T b) { if (b < a) { a = b; return 1; } return 0; } template bool chmax(T& a, const T b) { if (a < b) { a = b; return 1; } return 0; } template pair divmod(T a, T b) { assert(a >= 0 && b > 0); return make_pair(T(a / b), T(a % b)); } template T mod(T x, T m) { assert(m != 0); return T((x % m + m) % m); } template T ceil(T a, T b) { if (b < 0) { return ceil(T(-a), T(-b)); } assert(b > 0); return (a < 0 ? T(a / b) : T((a + b - 1) / b)); } template T floor(T a, T b) { if (b < 0) { return floor(T(-a), T(-b)); } assert(b > 0); return (a > 0 ? T(a / b) : T((a - b + 1) / b)); } template T powint(T x, T n) { assert(n >= 0); if (n == 0) { return T(1); }; T res = powint(x, T(n>>1)); res *= res; if (n & 1) { res *= x; } return res; } ll bitlen(ll b) { if (b <= 0) { return 0; } return (64LL - __builtin_clzll(b)); } ll digit_len(ll n) { assert(n >= 0); if (n == 0) { return 1; } ll sum = 0; while (n > 0) { sum++; n /= 10; } return sum; } ll digit_sum(ll n) { assert(n >= 0); ll sum = 0; while (n > 0) { sum += n % 10; n /= 10; } return sum; } ll digit_prod(ll n) { assert(n >= 0); if (n == 0) { return 0; } ll prod = 1; while (n > 0) { prod *= n % 10; n /= 10; } return prod; } ll xor_sum(ll x) { assert(0 <= x); switch (x % 4) { case 0: return x; case 1: return 1; case 2: return x ^ 1; case 3: return 0; } assert(false); } string toupper(const string& S) { string T(S); for (int i = 0; i < (int)T.size(); i++) { T[i] = toupper(T[i]); } return T; } string tolower(const string& S) { string T(S); for (int i = 0; i < (int)T.size(); i++) { T[i] = tolower(T[i]); } return T; } int a2i(const char& c) { assert(islower(c)); return (c - 'a'); } int A2i(const char& c) { assert(isupper(c)); return (c - 'A'); } int d2i(const char& d) { assert(isdigit(d)); return (d - '0'); } char i2a(const int& i) { assert(0 <= i && i < 26); return ('a' + i); } char i2A(const int& i) { assert(0 <= i && i < 26); return ('A' + i); } char i2d(const int& i) { assert(0 <= i && i <= 9); return ('0' + i); } using P = pair; using VP = vector

; using VVP = vector; using VS = vector; using VVS = vector; using VI = vector; using VVI = vector; using VVVI = vector; using VLL = vector; using VVLL = vector; using VVVLL = vector; using VB = vector; using VVB = vector; using VVVB = vector; using VD = vector; using VVD = vector; using VVVD = vector; using VLD = vector; using VVLD = vector; using VVVLD = vector; const ld EPS = 1e-10; const ld PI = acosl(-1.0); constexpr ll MOD = 1000000007; // constexpr ll MOD = 998244353; constexpr int inf = (1 << 30) - 1; // 1073741824 - 1 constexpr ll INF = (1LL << 62) - 1; // 4611686018427387904 - 1 // -------------------------------------------------------- // #include // using namespace atcoder; // References: // 『プログラミングコンテスト攻略のためのアルゴリズムとデータ構造』 // // // 微小値 EPS の誤差を許容して a, b が等価であるか判定する inline bool eq(ld a, ld b) { return fabs(a - b) < EPS; } // 点 (ベクトル) struct Point { ld x, y; Point(ld x = 0, ld y = 0): x(x), y(y) {} Point& operator+=(const Point& v) noexcept { x += v.x; y += v.y; return *this; } Point& operator-=(const Point& v) noexcept { x -= v.x; y -= v.y; return *this; } Point& operator*=(ld k) noexcept { x *= k; y *= k; return *this; } Point& operator/=(ld k) noexcept { x /= k; y /= k; return *this; } Point operator+(const Point& v) const noexcept { return Point(*this) += v; } Point operator-(const Point& v) const noexcept { return Point(*this) -= v; } Point operator*(ld k) const noexcept { return Point(*this) *= k; } Point operator/(ld k) const noexcept { return Point(*this) /= k; } ld norm() const noexcept { return x*x + y*y; } // ベクトルの大きさ ld abs() const noexcept { return sqrt(norm()); } // 原点からの距離,ベクトルの長さ // pair の要領で大小比較(x昇順 --> y昇順) bool operator < (const Point& p) const noexcept { return x != p.x ? x < p.x : y < p.y; } bool operator == (const Point& p) const noexcept { return eq(x, p.x) && eq(y, p.y); } }; // 直線 struct Line { Point p1, p2; Line(Point p1 = Point(), Point p2 = Point()): p1(p1), p2(p2) {} }; // 線分 struct Segment { Point p1, p2; Segment(Point p1 = Point(), Point p2 = Point()): p1(p1), p2(p2) {} }; // 円 struct Circle { Point c; // 中心 ld r; // 半径 Circle(Point c = Point(), ld r = 0.0): c(c), r(r) {} }; // 三角形 struct Triangle { Point p1, p2, p3; Triangle(Point p1 = Point(), Point p2 = Point(), Point p3 = Point()) : p1(p1), p2(p2), p3(p3) {} }; using Points = vector; using Polygon = vector; // 多角形 using Segments = vector; using Lines = vector; using Circles = vector; using Triangles = vector; ld norm(const Point& a) { return a.x*a.x + a.y*a.y; } // ベクトルの大きさ ld abs(const Point& a) { return sqrt(norm(a)); } // 原点からの距離,ベクトルの長さ ld dot(const Point& a, const Point& b) { return a.x*b.x + a.y*b.y; } // 内積 ld cross(const Point& a, const Point& b) { return a.x*b.y - a.y*b.x; } // 外積 // 射影: 直線 l に点 p から垂線を引いたときの交点を求める Point projection(const Line& l, const Point& p) { Point base = l.p2 - l.p1; ld t = dot(p - l.p1, base) / norm(base); return l.p1 + base * t; } // 射影: 線分 s に点 p から垂線を引いたときの交点を求める // ※ 交点が線分上に存在するかは別途判定する必要あり --> on_segment() Point projection(const Segment& s, const Point& p) { Point base = s.p2 - s.p1; ld t = dot(p - s.p1, base) / norm(base); return s.p1 + base * t; } // 反射: 直線 l を対称軸として点 p と線対称にある点を求める Point reflection(const Line& l, const Point& p) { return p + (projection(l, p) - p) * 2.0; } // 反射: 線分 s を対称軸として点 p と線対称にある点を求める // ※ 線分が交差するかは別途判定する必要あり --> is_intersected() Point reflection(const Segment& s, const Point& p) { return p + (projection(s, p) - p) * 2.0; } // 直交判定 (直線) bool is_orthogonal(const Line& a, const Line& b) { return eq(dot(a.p1 - a.p2, b.p1 - b.p2), 0.0); } // 直交判定 (線分) // ※ 線分の交差判定を先に行う必要あり --> is_intersected() bool is_orthogonal(const Segment& a, const Segment& b) { return eq(dot(a.p1 - a.p2, b.p1 - b.p2), 0.0); } // 平行判定 (直線) bool is_parallel(const Line& a, const Line& b) { return eq(cross(a.p1 - a.p2, b.p1 - b.p2), 0.0); } // 平行判定 (線分) bool is_parallel(const Segment& a, const Segment& b) { return eq(cross(a.p1 - a.p2, b.p1 - b.p2), 0.0); } constexpr int COUNTER_CLOCKWISE = 1; // 反時計回り (左回り) constexpr int CLOCKWISE = -1; // 時計回り (右回り) constexpr int ONLINE_BACK = 2; // 線分後方 constexpr int ONLINE_FRONT = -2; // 線分前方 constexpr int ON_SEGMENT = 0; // 線分上 // 3点 p0,p1,p2 の位置関係を線分 p1 - p0 を基準にして求める int ccw(const Point& p0, const Point& p1, const Point& p2) { Point a = p1 - p0; Point b = p2 - p0; if (cross(a, b) > EPS) return COUNTER_CLOCKWISE; // p0 -> p1, 反時計回りの方向に p2 if (cross(a, b) < -EPS) return CLOCKWISE; // p0 -> p1, 時計回りの方向に p2 if (dot(a, b) < -EPS) return ONLINE_BACK; // p2 -> p0 -> p1 の順で直線上に p2 if (a.norm() < b.norm()) return ONLINE_FRONT; // p0 -> p1 -> p2 の順で直線上に p2 return ON_SEGMENT; // p0 -> p2 -> p1 の順で線分上に p2 } // 線分の交差判定 // 端点だけ重なる場合や平行に重なる場合も交差しているとみなす bool is_intersected(const Point& p1, const Point& p2, const Point& p3, const Point& p4) { return (ccw(p1,p2,p3) * ccw(p1,p2,p4) <= 0) && (ccw(p3,p4,p1) * ccw(p3,p4,p2) <= 0); } // 線分の交差判定 // 端点だけ重なる場合や平行に重なる場合も交差しているとみなす bool is_intersected(const Segment& s1, const Segment& s2) { return is_intersected(s1.p1, s1.p2, s2.p1, s2.p2); } // 直線と線分の交差判定 // 端点だけ重なる場合や平行に重なる場合も交差しているとみなす bool is_intersected(const Line& l, const Segment& s) { return (ccw(l.p1, l.p2, s.p1) * ccw(l.p1, l.p2, s.p2) <= 0); } // 線分上に点が存在するか判定 bool on_segment(const Segment& s, const Point& p) { return ccw(s.p1, s.p2, p) == ON_SEGMENT; } // 線分の交点を求める // 線分 s1,s2 が交点を持ち平行に重なっていないことを想定 Point cross_point(const Segment& s1, const Segment& s2) { Point base = s2.p2 - s2.p1; ld d1 = abs(cross(base, s1.p1 - s2.p1)); ld d2 = abs(cross(base, s1.p2 - s2.p1)); assert(EPS < d1 || EPS < d2); return s1.p1 + (s1.p2 - s1.p1) * (d1 / (d1 + d2)); } // 直線の交点を求める // 直線 l1,l2 が交点を持ち平行に重なっていないことを想定 Point cross_point(const Line& l1, const Line& l2) { Point a = l1.p1, b = l1.p2; Point c = l2.p1, d = l2.p2; return a + (b - a) * cross(c - a, d - c) / cross(b - a, d - c); } // 直線と線分の交点を求める // 直線 l と線分 s が交点を持ち平行に重なっていないことを想定 Point cross_point(const Line& l, const Segment& s) { Point a = l.p1, b = l.p2; Point c = s.p1, d = s.p2; return a + (b - a) * cross(c - a, d - c) / cross(b - a, d - c); } // 2点間の距離 ld distance(const Point& a, const Point& b) { return abs(a - b); } // 点と直線の距離 ld distance(const Point& p, const Line& l) { return abs(p - projection(l, p)); } // 点と線分の距離 ld distance(const Point& p, const Segment& s) { Point m = projection(s, p); if (on_segment(s, m)) return abs(p - m); return min(abs(s.p1 - p), abs(s.p2 - p)); } // 線分と線分の距離 ld distance(const Segment& a, const Segment& b) { if (is_intersected(a, b)) return 0.0; return min({distance(a.p1, b), distance(a.p2, b), distance(b.p1, a), distance(b.p2, a)}); } // 円周上に点が存在するか判定 bool on_circle(const Circle& c, const Point& p) { return eq(c.r, distance(c.c, p)); } // 円内部に点が存在するか判定 // 円周上に点が存在する場合は false を返す bool in_circle(const Circle& c, const Point& p) { if (on_circle(c, p)) return false; return distance(c.c, p) < c.r; } // 円内部に線分が存在するか判定 // 円周上に端点が存在する場合は false を返す bool in_circle(const Circle& c, const Segment& s) { return in_circle(c, s.p1) && in_circle(c, s.p2); } // 円と直線が接するか判定 bool on_circle(const Circle& c, const Line& l) { return eq(c.r, distance(c.c, l)); } // 直線上に点が存在するか判定 bool on_line(const Line& l, const Point& p) { return eq(distance(p, l), 0.0); } // 円における点の内包判定 // 0: 外側, 1: 円周上, 2: 内側 int point_containment(const Circle& c, const Point& p) { if (in_circle(c, p)) return 2; if (on_circle(c, p)) return 1; return 0; } // 円と直線の交差判定 // 接している場合も交差しているとみなす bool is_intersected(const Circle& c, const Line& l) { if (on_circle(c, l)) return true; return distance(c.c, l) < c.r; } // 円の交差判定 // 2つの円の共通接線の数を計算する // 4本: 離れている // 3本: 外接 // 2本: 2点交差 // 1本: 内接 // 0本: 内包 int circle_intersection(const Circle& c1, const Circle& c2) { int n; // 共通接線の数 ld d = distance(c1.c, c2.c); ld r1 = c1.r, r2 = c2.r; if (r1 > r2) swap(r1, r2); // r1 <= r2 if (r1 + r2 < d) { n = 4; } else if (eq(r1 + r2, d)) { n = 3; } else if (d + r1 > r2) { n = 2; } else if (eq(d + r1, r2)) { n = 1; } else { // d + r1 < r2 n = 0; } return n; } // 円の交差判定 // 外接・2点交差・内接の場合に交差していると判定する bool is_intersected(const Circle& c1, const Circle& c2) { int n = circle_intersection(c1, c2); return (n == 3 || n == 2 || n == 1) ? true : false; } // 円と直線の交点を求める // 円と直線が交差していることを想定 // ただし,接している場合は接点が2つ返される pair cross_point(const Circle& c, const Line& l) { assert(is_intersected(c, l)); Point pr = projection(l, c.c); Point e = (l.p2 - l.p1) / abs(l.p2 - l.p1); ld base = sqrt(c.r * c.r - norm(pr - c.c)); return make_pair(pr + e * base, pr - e * base); } // 偏角を求める (argument) ld arg(const Point& p) { return atan2(p.y, p.x); } // 極座標系から直交座標系に変換する Point polar2carte(const ld& r, const ld& theta) { return Point(r * cos(theta), r * sin(theta)); } // 円と円の交点を求める // 円と円が交差していることを想定 // ただし,接している場合は接点が2つ返される pair cross_point(const Circle& c1, const Circle& c2) { assert(is_intersected(c1, c2)); ld d = distance(c1.c, c2.c); ld a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d)); // 余弦定理 ld t = arg(c2.c - c1.c); return make_pair(c1.c + polar2carte(c1.r, t + a), c1.c + polar2carte(c1.r, t - a)); } // 点 p を通る円 c の接線を求める // 点 p が円 c の外側に存在することを想定 pair tangent_point(const Point& p, const Circle& c) { // 三平方の定理より点 p と接点の距離 d が得られるため, // 点 p を中心とする半径 d の円と円 c の交点を求めればよい ld d = sqrt(norm(c.c - p) - c.r * c.r); return cross_point(c, Circle(p, d)); } // 2つの円の共通接線を求める // - (c1 の接点, c2 の接点) の順で Line を列挙する // Reference: // // // → 共通内接線/共通外接線の数式が逆になっている(コードは合ってる) Lines common_tangent(const Circle& c1, const Circle& c2) { ld d = distance(c1.c, c2.c); ld x1 = c1.c.x, y1 = c1.c.y, r1 = c1.r; ld x2 = c2.c.x, y2 = c2.c.y, r2 = c2.r; ld xd = x2 - x1, yd = y2 - y1; ld k1 = r1 / d, k2 = r2 / d; Lines lines; // 共通外接線 if (abs(r1 - r2) <= d) { ld c = (r1 - r2) / d; // cos ld s = sqrt(1 - c*c); // sin if (eq(abs(r1 - r2), d)) { assert(eq(s, 0)); Point p1(x1 + k1 * c*xd, y1 + k1 * c*yd); Point p2(x2 + k2 * c*xd, y2 + k2 * c*yd); lines.emplace_back(p1, p2); } else { Point p1a(x1 + k1 * (c*xd - s*yd), y1 + k1 * ( s*xd + c*yd)); Point p1b(x1 + k1 * (c*xd + s*yd), y1 + k1 * (-s*xd + c*yd)); Point p2a(x2 + k2 * (c*xd - s*yd), y2 + k2 * ( s*xd + c*yd)); Point p2b(x2 + k2 * (c*xd + s*yd), y2 + k2 * (-s*xd + c*yd)); lines.emplace_back(p1a, p2a); lines.emplace_back(p1b, p2b); } } // 共通内接線 if (r1 + r2 <= d) { ld c = (r1 + r2) / d; // cos ld s = sqrt(1 - c*c); // sin if (eq(r1 + r2, d)) { assert(eq(s, 0)); Point p1(x1 + k1 * c*xd, y1 + k1 * c*yd); Point p2(x2 + k2 * c*xd, y2 + k2 * c*yd); lines.emplace_back(p1, p2); } else { Point p1a(x1 + k1 * (c*xd - s*yd), y1 + k1 * ( s*xd + c*yd)); Point p1b(x1 + k1 * (c*xd + s*yd), y1 + k1 * (-s*xd + c*yd)); Point p2a(x2 - k2 * (c*xd - s*yd), y2 - k2 * ( s*xd + c*yd)); Point p2b(x2 - k2 * (c*xd + s*yd), y2 - k2 * (-s*xd + c*yd)); lines.emplace_back(p1a, p2a); lines.emplace_back(p1b, p2b); } } return lines; } // 円と多角形の共通部分を求める // TODO // 2つの円の共通部分の面積を求める // Reference: https://tjkendev.github.io/procon-library/python/geometry/circles_intersection_area.html ld area_of_intersection_of_two_circles(const Circle& c1, const Circle& c2) { ld d = distance(c1.c, c2.c); ld dd = d*d; ld r1 = c1.r; ld r2 = c2.r; int n_common_tangent = circle_intersection(c1, c2); // 離れている or 外接 if (3 <= n_common_tangent) return 0; // 内包 or 内接 if (n_common_tangent <= 1) { ld r = min(r1, r2); return PI*r*r; } ld p1 = r1*r1 - r2*r2 + dd; ld p2 = r2*r2 - r1*r1 + dd; ld S1 = r1*r1 * atan2(sqrt(dd*r1*r1*4 - p1*p1), p1); ld S2 = r2*r2 * atan2(sqrt(dd*r2*r2*4 - p2*p2), p2); ld S0 = sqrt(dd*r1*r1*4 - p1*p1) / 2; return -S0 + S1 + S2; } // 2点を直径とする円を求める Circle circle_with_2pt_as_diameter(const Point& p1, const Point& p2) { return Circle((p1 + p2) / 2, distance(p1, p2) / 2); } // 三角形の面積を求める ld area_of_triangle(const Point& A, const Point& B, const Point& C) { ld a = abs(B - C); ld b = abs(C - A); ld c = abs(A - B); ld s = (a + b + c) / 2; return sqrt(s*(s-a)*(s-b)*(s-c)); } // 三角形の面積を求める ld area_of_triangle(const Triangle& t) { return area_of_triangle(t.p1, t.p2, t.p3); } // 3 点が同一直線上にあるか判定する bool are_on_same_line(const Point& p1, const Point& p2, const Point& p3) { return abs(ccw(p1, p2, p3)) != 1; } // 三角形の内接円を求める // - 3点が同じ直線上に存在しない想定 // Circle incircle_of_triangle(const Point& A, const Point& B, const Point& C) { assert(not are_on_same_line(A, B, C)); // 3点が同じ直線上に存在しない確認 ld a = abs(B - C); ld b = abs(C - A); ld c = abs(A - B); ld d = a + b + c; ld s = d / 2; ld S = sqrt(s*(s-a)*(s-b)*(s-c)); ld r = 2*S / d; Point I = A*(a/d) + B*(b/d) + C*(c/d); return Circle(I, r); } // 三角形の外接円を求める // - 3点が同じ直線上に存在しない想定 // Circle circumscribed_circle_of_triangle(const Point& A, const Point& B, const Point& C) { assert(not are_on_same_line(A, B, C)); // 3点が同じ直線上に存在しない確認 ld a2 = norm(B - C); ld b2 = norm(C - A); ld c2 = norm(A - B); ld x = a2 * (b2 + c2 - a2); ld y = b2 * (c2 + a2 - b2); ld z = c2 * (a2 + b2 - c2); ld w = x + y + z; // 情報落ち対策のため x/w などとしている Point O = A*(x/w) + B*(y/w) + C*(z/w); // 外心 ld R = distance(O, A); // 外接円の半径 return Circle(O, R); } // 多角形の面積を求める // - 3 ≦ 頂点数 // - 頂点列が反時計回りであることを想定 // - 多角形が自己交差していないことを想定 // ld area_of_polygon(const Polygon& P) { int N = P.size(); assert(3 <= N); ld area = 0.0; for (int i = 0; i < N; i++) { area += cross(P[i], P[(i+1) % N]); // N番目の次は1番目 } return area / 2; } // 凸性判定 (多角形が凸多角形であるか判定) // - 3 ≦ 頂点数 // - 頂点列が反時計回りであることを想定 // - 多角形が自己交差していないことを想定 bool is_convex(const Polygon& P) { int N = P.size(); assert(3 <= N); for (int i = 0; i < N; i++) { if (ccw(P[(i-1+N) % N], P[i], P[(i+1+N) % N]) == CLOCKWISE) return false; } return true; } // 多角形における点の内包判定 // - 0: 外側, 1: 線分上, 2: 内側 int point_containment(const Polygon& g, const Point& p) { int N = g.size(); bool x = false; for (int i = 0; i < N; i++) { Point a = g[i] - p; Point b = g[(i+1) % N] - p; // N番目の次は1番目 if (abs(cross(a, b)) < EPS && dot(a, b) < EPS) return 1; // 線分上 if (a.y > b.y) swap(a, b); if (a.y < EPS && EPS < b.y && cross(a, b) > EPS) x = !x; // 半直線との交差回数の偶奇 } return (x ? 2 : 0); } // 凸包を求める (Andrew's Monotone Chain) // - 下記が満たされていることを想定 // - 3 ≦ 頂点数 // - 全点が一つの直線上に存在しない // -> 適当に2点取って直線を求め、残りの全点がその直線上に存在すれば NG // - on_edge: 凸包の辺上の点を含めるか // - 凸包の頂点の順序は最も左の頂点から反時計回り (左回り) Polygon convex_hull(Polygon P, bool on_edge = true) { int N = P.size(); assert(3 <= N); sort(P.begin(), P.end()); // x昇順 --> y昇順 Polygon ch(2*N,{-1,-1}); int k = 0; if (on_edge) { // 上包 (upper hull) for (int i = 0; i < N; ch[k++] = P[i++]) { while (k >= 2 && ccw(ch[k-2], ch[k-1], P[i]) == COUNTER_CLOCKWISE) k--; } // 下包 (lower hull) const int t = k + 1; for (int i = N-2; 0 <= i; ch[k++] = P[i--]) { while (k >= t && ccw(ch[k-2], ch[k-1], P[i]) == COUNTER_CLOCKWISE) k--; } } else { // 上包 (upper hull) for (int i = 0; i < N; ch[k++] = P[i++]) { while (k >= 2 && ccw(ch[k-2], ch[k-1], P[i]) != CLOCKWISE) k--; } // 下包 (lower hull) const int t = k + 1; for (int i = N-2; 0 <= i; ch[k++] = P[i--]) { while (k >= t && ccw(ch[k-2], ch[k-1], P[i]) != CLOCKWISE) k--; } } ch.resize(k-1); /** TODO: 最初から反時計回りで計算 **/ reverse(ch.begin() + 1, ch.end()); // 反時計回りに変換 return ch; } // 凸多角形を直線で切断したとき左側にできる凸多角形を求める // - 3 ≦ 頂点数 // - 頂点列が反時計回りであることを想定 // - 直線上に位置する頂点も結果に含まれる // - 求めた凸多角形の頂点数が 0〜2 になる可能性あり (面積計算等で注意) Polygon convex_cut(const Polygon& P, const Line& l) { int N = P.size(); assert(3 <= N); Polygon cc; Point p1 = l.p1; Point p2 = l.p2; for (int i = 0; i < N; i++) { Point a = P[i]; Point b = P[(i+1) % N]; // N番目の次は1番目 Segment s(a, b); // 1. 点 a が直線の左側に位置すれば cc に追加 (直線上にある場合も含む) if (ccw(p1, p2, a) != CLOCKWISE) { cc.push_back(a); } // 2. 線分 ab が直線と交差すれば交点を cc に追加 (端点の場合は含めない) if (!on_line(l, a) && !on_line(l, b) && is_intersected(l, s)) { cc.push_back(cross_point(l, Line(a, b))); } } return cc; } // 凸多角形の直径を求める(最遠点対間距離) // - Q が凸多角形であることを想定 (is_convex(Q) が true) // -> convex_hull(P, true) で計算して引数に指定 ld convex_diameter(Polygon Q) { int N = Q.size(); if (N == 2) { return distance(Q[0], Q[1]); } int i = 0, j = 0; // 最遠点対 for (int k = 1; k < N; k++) { // x軸方向の最遠点対で初期化 if (Q[k] < Q[i]) i = k; // 最も左側 if (Q[j] < Q[k]) j = k; // 最も右側 } ld res = 0; int si = i, sj = j; // 終了判定用 while (i != sj || j != si) { // 反時計回りに180度まで回転しながら解析 res = max(res, distance(Q[i], Q[j])); // 符号付き面積を利用して次の最遠点対が Q_{i+1}, Q_{j+1} のどちらか決定 if (cross(Q[(i+1) % N] - Q[i], Q[(j+1) % N] - Q[j]) < 0) { i = (i + 1) % N; } else { j = (j + 1) % N; } } return res; } // 最近点対間距離を求める ld closest_pair(Points P) { int N = P.size(); assert(2 <= N); sort(P.begin(), P.end()); auto compare_y = [](const Point& a, const Point& b) -> bool { return a.y < b.y; }; constexpr ld DINF = 1e18; Points B(N); // x座標を左右に分ける直線付近の点集合を保管用に使い回す配列 auto rec = [&](auto&& self, int l, int r) -> ld { if (r - l <= 1) return DINF; int m = (l + r) / 2; // x座標を左右に分ける直線のインデックス ld x = P[m].x; // x座標を左右に分ける直線のx座標 ld d = min(self(self, l, m), self(self, m, r)); // 最近点対間距離 inplace_merge(P.begin() + l, P.begin() + m, P.begin() + r, compare_y); int k = 0; for (int i = l; i < r; i++) { if (abs(P[i].x - x) >= d) continue; for (int j = k - 1; 0 <= j; j--) { ld dx = P[i].x - B[j].x; ld dy = P[i].y - B[j].y; if (dy >= d) break; d = min(d, sqrt(dx*dx + dy*dy)); } B[k++] = P[i]; } return d; }; return rec(rec, 0, N); } // 多角形の頂点列が反時計回りか判定する // - 符号付き面積の正負で判定する bool is_counter_clockwise(const Polygon& P) { return area_of_polygon(P) > 0.0; } // 三角形の頂点列が反時計回りか判定する // - 符号付き面積の正負で判定する bool is_counter_clockwise(const Triangle& t) { return area_of_triangle(t) > 0.0; } // 三角形における点の内包判定 // - 0: 外側, 1: 線分上, 2: 内側 int point_containment(const Triangle& t, const Point& p) { int c12 = ccw(t.p1, t.p2, p); int c23 = ccw(t.p2, t.p3, p); int c31 = ccw(t.p3, t.p1, p); if (c12 == 0 || c23 == 0 || c31 == 0) { return 1; } // 辺上 if (c12 == +1 && c23 == +1 && c31 == +1) { return 2; } // 時計回り if (c12 == -1 && c23 == -1 && c31 == -1) { return 2; } // 反時計回り return 0; } // 多角形 P において三角形 (P[l],P[m],P[r]) が耳であるか判定する bool is_ear(const Polygon& P, int l, int m, int r) { Triangle t(P[l], P[m], P[r]); // 明らかに ear でないケースを除外(この 3 点が mouth or 一直線) if (ccw(t.p1, t.p2, t.p3) <= 0) { return false; } int N = P.size(); for (int i = 0; i < N; i++) { if (i == l || i == m || i == r) { continue; } if (point_containment(t, P[i])) { continue; } } return true; } // 単純多角形を三角形分割する(耳分解) // - 3 ≦ 頂点数 // - 頂点列が反時計回りであることを想定 // - O(N^2) : ear clipping method によるナイーブな実装 // - Reference: https://web.archive.org/web/20200222054711/http://www.prefield.com/algorithm/geometry/triangulate.html Triangles triangulation_of_polygon(const Polygon& P) { int N = P.size(); assert(2 <= N); vector L(N), R(N); // 左隣・右隣 for (int i = 0; i < N; i++) { L[i] = (i - 1 + N) % N; R[i] = (i + 1 + N) % N; } Triangles ts; int m = 0; // いま見ている3点の真ん中 while ((int)ts.size() < N-2) { m = R[m]; while (is_ear(P, L[m], m, R[m])) { ts.emplace_back(P[L[m]], P[m], P[R[m]]); // 真ん中の m を闇に葬る R[L[m]] = R[m]; L[R[m]] = L[m]; } } return ts; } // 三角形と直線の共通部分の長さを求める ld length_of_intersection_of_triangle_and_line(const Triangle& t, const Line& l) { Segments ss; ss.emplace_back(t.p1, t.p2); ss.emplace_back(t.p2, t.p3); ss.emplace_back(t.p3, t.p1); Points ps; for (const auto& s : ss) { // 直線上に辺が存在 if (on_segment(s, l.p1) && on_segment(s, l.p2)) { return distance(s.p1, s.p2); } // 線分と直線の交点 if (is_intersected(l, s)) { ps.emplace_back(cross_point(l, s)); } } // 交点は高々 2 種類しかないため全探索 ld res = 0; int n = ps.size(); for (int i = 0; i < n-1; i++) { for (int j = i+1; j < n; j++) { res = max(res, distance(ps[i], ps[j])); } } return res; } // References: // // // Disjoint-set data structure (Union Find) // struct dsu { public: dsu() : N(0) {} explicit dsu(int n) : N(n), parent_or_size(n, -1), n_edge(n, 0) {} // 辺 (a, b) を張ってマージ成否を返す : amortized O(α(N)) bool merge(int a, int b) { assert(0 <= a && a < N); assert(0 <= b && b < N); int x = leader(a), y = leader(b); if (x == y) { n_edge[x]++; return false; } if (-parent_or_size[x] < -parent_or_size[y]) { swap(x, y); } parent_or_size[x] += parent_or_size[y]; parent_or_size[y] = x; n_edge[x] += n_edge[y] + 1; return true; } // 頂点 a, b が連結か判定する : amortized O(α(N)) bool same(int a, int b) { assert(0 <= a && a < N); assert(0 <= b && b < N); return leader(a) == leader(b); } // 頂点 a の属する連結成分のルートを返す : amortized O(α(N)) int leader(int a) { assert(0 <= a && a < N); if (parent_or_size[a] < 0) { return a; } return parent_or_size[a] = leader(parent_or_size[a]); } // 頂点 a が属する連結成分のサイズを返す : amortized O(α(N)) int size(int a) { assert(0 <= a && a < N); return -parent_or_size[leader(a)]; } // a が属する連結成分の辺の数を返す : amortized O(α(N)) int size_e(int a) { assert(0 <= a && a < N); return n_edge[leader(a)]; } // 「一つの連結成分の頂点番号リスト」のリストを返す : O(N) vector> groups() { vector leader_buf(N), group_size(N); for (int i = 0; i < N; i++) { leader_buf[i] = leader(i); group_size[leader_buf[i]]++; } vector> result(N); for (int i = 0; i < N; i++) { result[i].reserve(group_size[i]); } for (int i = 0; i < N; i++) { result[leader_buf[i]].push_back(i); } result.erase( remove_if(result.begin(), result.end(), [&](const vector& v) { return v.empty(); }), result.end()); return result; } private: int N; // [x < 0] -x が連結成分のサイズに対応 // [0 <= x] x が parent に対応 vector parent_or_size; vector n_edge; }; int main() { ios::sync_with_stdio(false); cin.tie(nullptr); cout << fixed << setprecision(15); ll N; input(N); VLL X(N),Y(N); REP(i,N) input(X[i], Y[i]); map mp; REP(i,N) mp[{X[i],Y[i]}] = i; dsu uf(N); REP(i,N) { ll x1 = X[i], y1 = Y[i]; FOR(x,-10,10+1) REP(y,10+1) { if (x*x + y*y > 100) break; ll x2 = x1 + x, y2 = y1 + y; if (mp.count({x2, y2})) { ll j = mp[{x2, y2}]; uf.merge(i, j); } } } ld ans = 1; for (auto& g : uf.groups()) { int n = SZ(g); Polygon P(n); REP(i,n) P[i] = {ld(X[g[i]]), ld(Y[g[i]])}; if (n == 1) { chmax(ans, (ld)2.0); } else if (n == 2) { chmax(ans, distance(P[0], P[1]) + 2); } else { bool all_on_same_line = true; FOR(i,2,n) { all_on_same_line &= are_on_same_line(P[0], P[1], P[i]); } if (all_on_same_line) { SORT(P); ld d = distance(P[0], P[n-1]); chmax(ans, d + 2); } else { auto Q = convex_hull(P); auto d = convex_diameter(Q); chmax(ans, d + 2); } } } print(ans); return 0; }