/*
Catalin-Stefan Tiseanu
Pre-written code is assembled from various sources found online.
Big thanks to the community for that !
*/
// Pre-written code below
using namespace std;
#include <algorithm>
#include <bitset>
#include <cassert>
#include <cctype>
#include <climits>
#include <cmath>
#include <complex>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <deque>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <utility>
#include <vector>
#include <unordered_map>
//#define DEBUG
#ifndef NDEBUG
# define assert_msg(condition, message) \
do { \
if (! (condition)) { \
std::cerr << "Assertion `" #condition "` failed in " << __FILE__ \
<< " line " << __LINE__ << ": " << message << std::endl; \
std::exit(EXIT_FAILURE); \
} \
} while (false)
#else
# define ASSERT(condition, message) do { } while (false)
#endif
#ifdef DEBUG
#define debug(args...) {dbg,args; cerr<<endl;}
#else
#define debug(args...) // Just strip off all debug tokens
#endif
template<class T> std::ostream& operator <<(std::ostream& stream, const vector<T> & v) {
for (auto el : v)
stream << el << " ";
return stream;
}
template<class T> std::ostream& operator <<(std::ostream& stream, const vector<vector<T>> & v) {
for (auto line : v) {
for (auto el : line)
stream << el << " ";
stream << endl;
}
return stream;
}
template<class T, class U> std::ostream& operator <<(std::ostream& stream, const pair<U, T> & p) {
stream << "( " << p.first << ", " << p.second << " )";
return stream;
}
class debugger {
public:
template<class T> void output (const T& v) {
cerr << v << " ";
}
template<class T> debugger& operator , (const T& v) {
output(v);
return *this;
}
} dbg;
/*
template<> void debugger::output (const <vector<vector<int>> & v) {
for (auto line : v) {
for (auto el : line)
cerr << v << " ";
cerr << endl;
}
}
*/
////////////////
// templates //
////////////////
// general
template<typename T> int size(const T& c) { return int(c.size()); }
template<typename T> T abs(T x) { return x < 0 ? -x : x; }
template<typename T> T sqr(T x) { return x*x; }
template<typename T> bool remin(T& x, T y) { if (x <= y) return false; x = y; return true; }
template<typename T> bool remax(T& x, T y) { if (x >= y) return false; x = y; return true; }
template<class T> inline vector<pair<T,int> > factorize(T n)//NOTES:factorize(
{vector<pair<T,int> > R;for (T i=2;n>1;){if (n%i==0){int C=0;for (;n%i==0;C++,n/=i);R.push_back(make_pair(i,C));}
i++;if (i>n/i) i=n;}if (n>1) R.push_back(make_pair(n,1));return R;}
// strings
template<typename T> string toStr(T x) { stringstream ss; ss << x; return ss.str(); }
// misc
#define EPS 1e-10
// types
typedef long long int64;
typedef unsigned long long uint64;
// shortcuts
#define all(_xx) _xx.begin(), _xx.end()
#define pb push_back
#define vi vector<int>
#define vvi vector<vector<int>>
#define vd vector<double>
#define vpii vector<pair<int,int> >
#define vpdd vector<pair<double,double> >
#define pii pair<int,int>
#define pdd pair<double, double>
#define pll pair<int64, int64>
#define mp(XX, YY) make_pair(XX, YY)
//#define fi first
//#define se second
#define X first
#define Y second
#define ll long long
#define SS stringstream
// for loops
#define re(II, NN) for (int II(0), _NN(NN); (II) < (NN); ++(II))
#define fod(II, XX, YY) for (int II(XX), _YY(YY); (II) >= (_YY); --(II))
#define fo(II, XX, YY) for (int II(XX), _YY(YY); (II) <= (_YY); ++(II))
#define foreach(it,c) for(__typeof((c).begin()) it=(c).begin();it!=(c).end();++it)
// Useful hardware instructions
#define bitcount __builtin_popcount
#define gcd __gcd
// Useful all around
#define checkbit(n,b) ( (n >> b) & 1)
#define DREP(a) sort(all(a));a.erase(unique(all(a)),a.end())
#define INDEX(arr,ind) (lower_bound(all(arr),ind)-arr.begin())
#define PAUSE cerr << "Press any key to continue ..."; char xxxx; scanf("%c", &xxxx);
#define fill(xx,val) memset(xx, val, sizeof(xx))
/////////////////// MISOF'S GEOMETRY LIBRARY ///////////////////////////////////////////////////
/* Surely won't compile with COORD_TYPE integer:
* size, normalize, left_side_angle
*
* (If compiles), won't work properly with COORD_TYPE integer:
* poly_area, center_of_mass, intersect_*, rotate, distance_*, circumcircle_center
*/
template<class T> void checkmin(T & a, const T & b) { if (b < a) a = b; }
template<class T> void checkmax(T & a, const T & b) { if (b > a) a = b; }
#define COORD_TYPE double // data type used to store the coordinates: %d || %lld || %lf || %Lf
typedef complex<COORD_TYPE> point;
typedef vector<point> polygon; // if used to store a polygon, don't repeat the first vertex
#define EPSILON (1e-4) // epsilon used for computations involving doubles ; dec to 1e-9 for %Lf
#ifndef M_PI // UVa judge doesn't know M_PI
#define M_PI 3.14159265358979323846
#endif
#define ABS(x) ({ __typeof(x) _tmp = (x); if (_tmp<0) _tmp = -_tmp; _tmp; })
// safe comparison with 0: is_zero, is_negative, is_positive, signum {{{
bool is_negative(int x) { return x<0; } bool is_zero(int x) { return x==0; } bool is_positive(int x) { return x>0; }
bool is_negative(long long x) { return x<0; } bool is_zero(long long x) { return x==0; } bool is_positive(long long x) { return x>0; }
bool is_negative(double x) { return x < -EPSILON; } bool is_zero(double x) { return fabs(x) <= EPSILON; } bool is_positive(double x) { return x > EPSILON; }
bool is_negative(long double x) { return x < -EPSILON; } bool is_zero(long double x) { return fabsl(x) <= EPSILON; } bool is_positive(long double x) { return x > EPSILON; }
template<class T> int signum(const T &A) { if (is_zero(A)) return 0; if (is_negative(A)) return -1; return 1; }
// }}}
// safe equality test for points {{{
template<class T> bool are_equal(const complex<T> &A, const complex<T> &B) { return is_zero(real(B)-real(A)) && is_zero(imag(B)-imag(A)); }
// }}}
// cross-product, dot_product, square_size (=dot_product(A,A)) for 2D vectors {{{
template<class T> T square_size(const complex<T> &A) { return real(A) * real(A) + imag(A) * imag(A); }
template<class T> T dot_product(const complex<T> &A, const complex<T> &B) { return real(A)*real(B) + imag(A)*imag(B); }
template<class T> T cross_product(const complex<T> &A, const complex<T> &B) { return real(A) * imag(B) - real(B) * imag(A); }
// }}}
// size, normalize for 2D real vectors {{{
template<class T> T size(const complex<T> &A) { return sqrt(real(A) * real(A) + imag(A) * imag(A)); }
template<class T> void normalize(complex<T> &A) { T Asize = size(A); A *= (1/Asize); }
// }}}
// safe colinearity and orientation tests: colinear, clockwise, counterclockwise {{{
template<class T> bool colinear(const complex<T> &A, const complex<T> &B, const complex<T> &C) { return is_zero( cross_product( B-A, C-A )); }
template<class T> bool colinear(const complex<T> &B, const complex<T> &C) { return is_zero( cross_product( B, C )); }
template<class T> bool clockwise(const complex<T> &A, const complex<T> &B, const complex<T> &C) { return is_negative( cross_product( B-A, C-A )); }
template<class T> bool clockwise(const complex<T> &B, const complex<T> &C) { return is_negative( cross_product( B, C )); }
template<class T> bool counterclockwise(const complex<T> &A, const complex<T> &B, const complex<T> &C) { return is_positive( cross_product( B-A, C-A )); }
template<class T> bool counterclockwise(const complex<T> &B, const complex<T> &C) { return is_positive( cross_product( B, C )); }
// }}}
// safe comparison function according to [argument,size]: compare_arg {{{
template<class T> bool compare_arg(const complex<T> &A, const complex<T> &B) {
// [0,0] ide uplne na zaciatok
if (are_equal(B,point(0,0))) return 0;
if (are_equal(A,point(0,0))) return 1;
// chceme poradie: pod osou x zlava doprava, kladna poloos, nad osou sprava dolava, zaporna poloos
int sgnA = signum(imag(A));
int sgnB = signum(imag(B));
if (sgnA == 0) if (signum(real(A))<0) sgnA = 2;
if (sgnB == 0) if (signum(real(B))<0) sgnB = 2;
if (sgnA != sgnB) return (sgnA < sgnB);
// v ramci polroviny sortime podla clockwise
if (counterclockwise(A,B)) return 1;
if (clockwise(A,B)) return 0;
// a ak sa este nerozhodlo, podla vzdialenosti, blizsie skor
return (square_size(A) < square_size(B));
}
// }}}
// safe comparison functions acc. to [x,y] and [y,x]: compare_XY, compare_YX {{{
template<class T> bool compare_XY(const complex<T> &A, const complex<T> &B) { if (!is_zero(real(A)-real(B))) return (is_negative(real(A)-real(B))); return (is_negative(imag(A)-imag(B))); }
template<class T> bool compare_YX(const complex<T> &A, const complex<T> &B) { if (!is_zero(imag(A)-imag(B))) return (is_negative(imag(A)-imag(B))); return (is_negative(real(A)-real(B))); }
// }}}
#ifdef LESS_THAN_COMPARES_LEXICOGRAPHICALLY
// bool operator < such that (a+bi) < (c+di) <==> [a,b] < [c,d] {{{
template<class T> bool operator < (const complex<T> &A, const complex<T> &B) { if (real(A) != real(B)) return (real(A) < real(B)); return (imag(A) < imag(B)); }
template<class T> bool operator < (const T &A, const complex<T> &B) { return complex<T>(A) < B; }
template<class T> bool operator < (const complex<T> &A, const T &B) { return A < complex<T>(B); }
template<class T> bool operator > (const complex<T> &A, const complex<T> &B) { return B < A; }
template<class T> bool operator > (const T &A, const complex<T> &B) { return B < A; }
template<class T> bool operator > (const complex<T> &A, const T &B) { return B < A; }
// }}}
#endif
// polygon area: twice_signed_poly_area, poly_area {{{
template<class T> T twice_signed_poly_area(const vector< complex<T> > &V) { T res = 0; for (unsigned int i=0; i<V.size(); i++) res += cross_product( V[i], V[(i+1)%V.size()] ); return res; }
template<class T> T poly_area(const vector< complex<T> > &V) { return fabs(0.5 * twice_signed_poly_area(V)); }
// }}}
// compute the center of mass of a polygon: center_of_mass {{{
template<class T> complex<T> center_of_mass(const vector< complex<T> > &V ) {
complex<T> sum(0,0);
T twice_area = 0.0;
for (unsigned i=0; i<V.size(); i++) {
sum += cross_product( V[i], V[(i+1)%V.size()] ) * ( V[i] + V[(i+1)%V.size()] );
twice_area += cross_product( V[i], V[(i+1)%V.size()] );
}
sum *= 1.0 / (3.0 * twice_area);
return sum;
}
// }}}
// compute a CCW convex hull with no unnecessary points: convex_hull {{{
template<class T> vector< complex<T> > convex_hull(vector< complex<T> > V) {
// handle boundary cases
if (V.size() == 2) if (are_equal(V[0],V[1])) V.pop_back();
if (V.size() <= 2) return V;
// find the bottommost point -- this can be optimized!
sort(V.begin(), V.end(), compare_YX<COORD_TYPE> );
complex<T> offset = V[0];
for (unsigned int i=0; i<V.size(); i++) V[i] -= offset;
sort(V.begin()+1, V.end(), compare_arg<COORD_TYPE> );
int count = 2;
vector<int> hull(V.size()+2);
for (unsigned int i=0; i<2; i++) hull[i]=i;
for (unsigned int i=2; i<V.size(); i++) {
while (count>=2 && !counterclockwise( V[hull[count-2]], V[hull[count-1]], V[i] ) ) count--;
hull[count++]=i;
}
vector< complex<T> > res;
for (int i=0; i<count; i++) res.push_back( V[hull[i]] + offset );
if (res.size()==2) if (are_equal(res[0],res[1])) res.pop_back();
return res;
}
// }}}
// safe test whether a point C \in [A,B], C \in (A,B): is_on_segment, is_inside_segment {{{
template<class T> bool is_on_segment(const complex<T> &A, const complex<T> &B, const complex<T> &C) {
if (!colinear(A,B,C)) return 0;
return ! is_positive( dot_product(A-C,B-C) );
}
template<class T> bool is_inside_segment(const complex<T> &A, const complex<T> &B, const complex<T> &C) {
if (!colinear(A,B,C)) return 0;
return is_negative( dot_product(A-C,B-C) );
}
// }}}
// winding number of a poly around a point (not on its boundary): winding_number(pt,poly) {{{
template<class T> int winding_number( const complex<T> &bod, const vector< complex<T> > &V ) {
int wn = 0;
for (unsigned int i=0; i<V.size(); i++) {
if (! is_positive( imag(V[i]) - imag(bod) )) {
if (is_positive( imag(V[(i+1)%V.size()]) - imag(bod) ))
if (counterclockwise( V[i], V[(i+1)%V.size()], bod ))
++wn;
} else {
if (! is_positive( imag(V[(i+1)%V.size()]) - imag(bod) ))
if (clockwise( V[i], V[(i+1)%V.size()], bod ))
--wn;
}
}
return wn;
}
// }}}
// test whether a point is inside a polygon: is_inside {{{
template<class T> int is_inside( const complex<T> &bod, const vector< complex<T> > &V ) {
for (unsigned int i=0; i<V.size(); i++) if (is_on_segment(V[i],V[(i+1)%V.size()],bod)) return true;
return winding_number(bod,V) != 0;
}
// }}}
// test whether a polygon is convex: is_poly_convex {{{
template<class T> bool is_poly_convex(const vector< complex<T> > &V) {
int cw=0, ccw=0;
int N = V.size();
for (unsigned i=0; i<V.size(); i++) if (clockwise(V[i],V[(i+1)%N],V[(i+2)%N])) { cw=1; break; }
for (unsigned i=0; i<V.size(); i++) if (counterclockwise(V[i],V[(i+1)%N],V[(i+2)%N])) { ccw=1; break; }
return (! (cw && ccw));
}
// }}}
// test whether a polygon is given in clockwise order: is_poly_clockwise {{{
template<class T> bool is_poly_clockwise(const vector< complex<T> > &V) {
return is_negative(twice_signed_poly_area(V));
}
// }}}
// intersect lines (A,B) and (C,D), return 2 pts if identical: intersect_line_line {{{
template<class T> vector< complex<T> > intersect_line_line(const complex<T> &A, const complex<T> &B, const complex<T> &C, const complex<T> &D) {
vector< complex<T> > res;
complex<T> U = B-A, V = D-C;
if (colinear(U,V)) { // identical or parallel
if (colinear(U,C-A)) { res.push_back(A); res.push_back(B); }
return res;
}
// one intersection point
T k = ( cross_product(C,V) - cross_product(A,V) ) / cross_product(U,V);
res.push_back(A + k*U);
return res;
}
// }}}
// intersect segments [A,B] and [C,D], may return endpoints of a segment: intersect_segment_segment {{{
template<class T> vector< complex<T> > intersect_segment_segment(const complex<T> &A, const complex<T> &B, const complex<T> &C, const complex<T> &D) {
vector< complex<T> > res;
complex<T> U = B-A, V = D-C, W = C-A, X = D-A;
if (colinear(U,V)) { // parallel
// check for degenerate cases
if (are_equal(A,B)) { if (is_on_segment(C,D,A)) res.push_back(A); return res; }
if (are_equal(C,D)) { if (is_on_segment(A,B,C)) res.push_back(C); return res; }
// two parallel segments
if (!colinear(U,W)) return res; // not on the same line
T ssU = square_size(U), ssW = square_size(W), ssX = square_size(X);
// does A coincide with C or D?
if (is_zero(ssW)) {
res.push_back(A);
if (is_negative(dot_product(U,X))) return res; // B, D on different sides
if (is_negative(ssU-ssX)) res.push_back(B); else res.push_back(D);
return res;
}
if (is_zero(ssX)) {
res.push_back(A);
if (is_negative(dot_product(U,W))) return res; // B, C on different sides
if (is_negative(ssU-ssW)) res.push_back(B); else res.push_back(C);
return res;
}
if (is_inside_segment(C,D,A)) {
// A is inside CD, intersection is [A,?]
res.push_back(A);
if (is_positive(dot_product(U,W))) { // B,C on the same side
if (is_negative(ssU-ssW)) res.push_back(B); else res.push_back(C);
return res;
} else { // B,D on the same side
if (is_negative(ssU-ssX)) res.push_back(B); else res.push_back(D);
return res;
}
} else {
// segment CD is strictly on one side of A
if (ssW < ssX) {
// A_C_D
if (is_negative(dot_product(U,W))) return res; // B is on the other side
if (is_negative(ssU-ssW)) return res; // B before C
if (are_equal(B,C)) { res.push_back(B); return res; }
res.push_back(C);
if (is_negative(ssU-ssX)) res.push_back(B); else res.push_back(D);
return res;
} else {
// A_D_C
if (is_negative(dot_product(U,X))) return res; // B is on the other side
if (is_negative(ssU-ssX)) return res; // B before D
if (are_equal(B,D)) { res.push_back(B); return res; }
res.push_back(D);
if (is_negative(ssU-ssW)) res.push_back(B); else res.push_back(C);
return res;
}
}
}
// not parallel, at most one intersection point
T k = ( cross_product(C,V) - cross_product(A,V) ) / cross_product(U,V);
complex<T> cand = A + k*U;
if (is_on_segment(A,B,cand) && is_on_segment(C,D,cand)) res.push_back(cand);
return res;
}
// }}}
// intersect a poly and a halfplane to the left of [A,B): intersect_poly_halfplane {{{
template<class T> vector< complex<T> > intersect_poly_halfplane(const vector< complex<T> > &V, const complex<T> &A, const complex<T> &B) {
int N = V.size();
vector< complex<T> > res;
if (N == 0) return res;
if (N == 1) if (! clockwise(A,B,V[0]) ) return V; else return res;
for (int i=0; i<N; i++) {
if (! clockwise(A,B,V[i])) res.push_back(V[i]);
int intersects = 0;
if (counterclockwise(A,B,V[i])) if (clockwise(A,B, V[(i+1)%N] )) intersects = 1;
if (clockwise(A,B,V[i])) if (counterclockwise(A,B, V[(i+1)%N] )) intersects = 1;
if (intersects) res.push_back( intersect_line_line(A,B,V[i], V[(i+1)%N] ));
}
return res;
}
// }}}
// intersect two convex polygons, result is CCW: intersect_cpoly_cpoly {{{
template<class T> vector< complex<T> > intersect_cpoly_cpoly(vector< complex<T> > V1, vector< complex<T> > V2) {
if (is_poly_clockwise(V1)) reverse(V1.begin(),V1.end());
if (is_poly_clockwise(V2)) reverse(V2.begin(),V2.end());
for (unsigned int i=0; i<V2.size(); i++) V1 = intersect_poly_halfplane(V1,V2[i],V2[(i+1) % V2.size() ]);
return V1;
}
// }}}
// intersect circles (C1,r1) and (C2,r2), 3 pts returned if \equiv: intersect_circle_circle {{{
template<class T> inline vector< complex<T> > intersect_circle_circle(const complex<T> &C1, T r1, const complex<T> &C2, T r2) {
vector< complex<T> > res;
if (are_equal(C1,C2)) {
// 2x the same point
if (is_zero(r1) && is_zero(r2)) { res.push_back(C1); return res; }
// identical circles -- return 3 points
if (is_zero(r1-r2)) { res.push_back(C1); res.push_back(C1); res.push_back(C1); return res; }
// no intersection
return res;
}
T d = sqrt(square_size(C2-C1));
// check for no intersection
if (is_positive(d-r1-r2) || is_positive(r1-r2-d) || is_positive(r2-r1-d)) return res;
// check for a single intersection
if (is_zero(d-r1-r2)) { res.push_back( (1.0/d) * (r1*C2 + r2*C1) ); return res; }
if (is_zero(r1-r2-d)) { res.push_back( C1 + (r1/d) * (C2-C1) ); return res; }
if (is_zero(r2-r1-d)) { res.push_back( C2 + (r2/d) * (C1-C2) ); return res; }
// general case: compute x and y offset of the intersections
T x = ( d*d - r2*r2 + r1*r1 ) / (2*d);
T y = sqrt( 4*d*d*r1*r1 - ( d*d - r2*r2 + r1*r1 )*( d*d - r2*r2 + r1*r1 ) ) / (2*d);
// I = (C1,C2) \cap chord, N = normal vector
complex<T> I = (1.0/d) * ( (d-x)*C1 + x*C2 );
complex<T> N( imag(C2-C1), -real(C2-C1) );
T Nsize = sqrt(square_size(N));
N = N * (1/Nsize);
// compute and return the points in lexicographic order
complex<T> I1 = I + y*N;
complex<T> I2 = I - y*N;
if (is_positive(real(I1)-real(I2))) swap(I1,I2);
if (is_zero(real(I1)-real(I2))) if (is_positive(imag(I1)-imag(I2))) swap(I1,I2);
res.push_back(I1);
res.push_back(I2);
return res;
}
// }}}
// compute distance: point A, line (B,C): distance_point_line {{{
template<class T> T distance_point_line(const complex<T> &A, const complex<T> &B, const complex<T> &C) {
complex<T> N ( imag(B-C), -real(B-C) );
normalize(N);
T tmp = dot_product(N,B-A);
return tmp;
}
// }}}
// compute distance: point A, segment [B,C]: distance_point_segment {{{
template<class T> T distance_point_segment(const complex<T> &A, const complex<T> &B, const complex<T> &C) {
T res = min( size(B-A), size(C-A) );
complex<T> N ( imag(B-C), -real(B-C) );
normalize(N);
T tmp = dot_product(N,B-A);
complex<T> paeta = A - tmp*N;
if (is_on_segment(B,C,paeta)) checkmin(res, fabs(tmp));
return res;
}
// }}}
// compute distance: line (A,B), line (C,D): distance_line_line {{{
template<class T> T distance_line_line(const complex<T> &A, const complex<T> &B, const complex<T> &C, const complex<T> &D) {
if (!colinear(B-A,D-C)) return 0.0;
return distance(A,C,D);
}
// }}}
// compute distance: segment [A,B], line (C,D): distance_segment_line {{{
template<class T> T distance_segment_line(const complex<T> &A, const complex<T> &B, const complex<T> &C, const complex<T> &D) {
if (! is_positive( cross_product(D-C,A-C) * cross_product(D-C,B-C) )) return 0; // they intersect
return min( distance(A,C,D), distance(B,C,D) );
}
// }}}
// compute distance: segment [A,B], segment [C,D]: distance_segment_segment
template<class T> T distance_segment_segment(const complex<T> &A, const complex<T> &B, const complex<T> &C, const complex<T> &D) {
vector< complex<T> > isect = intersect_segment_segment(A,B,C,D);
if (isect.size()) return 0; // they intersect
T res = min(min(distance_point_segment(A,C,D), distance_point_segment(B,C,D)),
min(distance_point_segment(C,A,B), distance_point_segment(D,A,B)));
return res;
}
// }}}
// circumcircle center for three points A,B,C (even colinear): {{{
template<class T> complex<T> circumcircle_center(const complex<T> &A, const complex<T> &B, const complex<T> &C) {
T a = 0, bx = 0, by = 0;
a += real(A) * imag(B) + real(B) * imag(C) + real(C) * imag(A);
a -= real(B) * imag(A) + real(C) * imag(B) + real(A) * imag(C);
if (is_zero(a)) { if (is_on_segment(A,B,C)) return 0.5*(A+B); if (is_on_segment(A,C,B)) return 0.5*(A+C); return 0.5*(B+C); }
bx += square_size(A) * imag(B) + square_size(B) * imag(C) + square_size(C) * imag(A);
bx -= square_size(B) * imag(A) + square_size(C) * imag(B) + square_size(A) * imag(C);
by -= square_size(A) * real(B) + square_size(B) * real(C) + square_size(C) * real(A);
by += square_size(B) * real(A) + square_size(C) * real(B) + square_size(A) * real(C);
return complex<T> ( bx / (2*a) , by / (2*a) );
}
// }}}
// rotate(point,center,CCW angle): {{{
template<class T> complex<T> rotate_point(const complex<T> &bod, const complex<T> &stred, T uhol) {
complex<T> mul(cos(uhol),sin(uhol)); return ((bod-stred)*mul)+stred;
}
// }}}
// rotate(poly,center,CCW angle): {{{
template<class T> vector< complex<T> > rotate_poly(vector< complex<T> > V, const complex<T> &stred, T uhol) {
for (unsigned int i=0; i<V.size(); i++) V[i] = rotate_point(V[i],stred,uhol); return V;
}
// }}}
// angle on the left side of B in polyline A->B->C: left_side_angle {{{
template<class T> T left_side_angle(const complex<T> &A, const complex<T> &B, const complex<T> &C) {
double a1 = atan2( imag(A-B), real(A-B) );
double a2 = atan2( imag(C-B), real(C-B) );
double u = a1 - a2;
if (u < 0) u += 2*M_PI;
if (u >= 2*M_PI) u -= 2*M_PI;
return u;
}
// }}}
// from here down added by Catalin-Stefan Tiseanu
// put a point X in Ax + By + C = 0
template<class T> T inline point_in_line(const complex<T> &P, const T &A, const T &B, const T &C) {
return A * real(P) + B * imag(P) + C;
}
// sign of a line
template<class T> inline int line_sign(const complex<T> & P, const T & A, const T & B, const T & C) {
return signum(A * real(P) + B * imag(P) + C);
}
// get coefficients of a line from two points (i.e Ax + By + C = 0)
template <class T> inline void get_line_coef(const complex<T> & P1, const complex<T> & P2, T &A, T &B, T &C ) {
A = imag(P1 - P2);
B = real(P2 - P1);
C = cross_product(P1, P2);
if (A < -EPSILON)
A = -A, B = -B, C = -C;
}
// get a point X on line (A, B) st. XB / XA = t
template<class T> complex<T> inline point_from_line(const complex<T> &A, const complex<T> &B, long double t) {
return A + (B - A) * t;
}
// get distance from point U to convex polygon ch
inline double distance_point_polygon(const point & U, const vector<point> & ch) {
double res = 1e20;
re (i, size(ch) - 1)
checkmin(res, distance_point_segment(U, ch[i], ch[i+1]));
checkmin(res, distance_point_segment(U, ch.front(), ch.back()));
/*
for (auto vertex : ch)
checkmin(res, size(vertex - U));
*/
return res;
}
//
/// END MISOF GEOMETRY LIBRARY
struct Scanner {
char* curPos;
const static int BUF_SIZE = 2*(1<<20);
char input[BUF_SIZE];
FILE * fin;
void init(FILE * _fin) {
fin = _fin;
fread(input, 1, sizeof(input), fin);
curPos = input;
}
Scanner() {;}
inline void ensureCapacity() {
int size = input + BUF_SIZE - curPos;
if (size < 100) {
memcpy(input, curPos, size);
fread(input + size, 1, BUF_SIZE - size, fin);
curPos = input;
}
}
inline int nextInt() {
ensureCapacity();
while (!isdigit(*curPos))
++curPos;
int res = 0;
while (*curPos >= '0' && *curPos <= '9')
res = res * 10 + (*(curPos++) & 15);
return res;
}
inline char nextChar() {
ensureCapacity();
while (*curPos <= ' ')
++curPos;
return *(curPos++);
}
} Sc;
const string filein = "oypara.in";
const string fileout = "oypara.out";
#define MAX_N 100111
int N;
polygon lower, upper, lower_hull, upper_hull;
// used to test if a line throgh U and V intersects the interior of the polygon ch
inline int same_side(const point & U, const point & V, const polygon & ch, int & side) {
int neg(0), pos(0);
side = 0;
double A, B, C;
get_line_coef(U, V, A, B, C);
re (i, size(ch)) {
int s = signum(line_sign(ch[i], A, B, C));
if (is_negative(s)) {++neg; side = -1; }
if (is_positive(s)) {++pos; side = 1; }
if (pos && neg)
return 0;
}
return 1;
}
inline int ternary_search(int l, int r, const polygon & from, const polygon & to) {
int o = (2 * l + r) / 3;
int p = (l + 2 * r) / 3;
//debug(l, r, o, p);
if (l + 3 >= r) {
int best_index = 1;
double res = 1e20;
fo (i, l, r) {
double cur_dist = distance_point_polygon(from[i], to);
if (cur_dist < res) {
res = cur_dist;
best_index = i;
}
}
return best_index;
}
if (distance_point_polygon(from[o], to) < distance_point_polygon(from[p], to))
return ternary_search(l, p, from, to);
else
return ternary_search(o, r, from, to);
}
// index of the closest point in <from> to the polygon <to>
/*
inline int closest_point_to_polygon(const polygon & from, const polygon & to) {
double res = 1e20;
int best_index = -1;
re (i, size(from)) {
double cur_dist = distance_point_polygon(from[i], to);
//debug("dist:", from[i], "to:", to, cur_dist);
//debug("dist", from[i], cur_dist);
if (cur_dist < res) {
res = cur_dist;
best_index = i;
}
}
return best_index;
}
*/
inline int closest_point_to_polygon(const polygon & from, const polygon & to) {
return ternary_search(0, size(from) - 1, from, to);
}
inline bool try_pair(const point & A, const point & B, const polygon & p1, const polygon &p2, pair<point, point> & line) {
debug("try ", A, B);
if (are_equal(A, B))
return false;
int s1(0), s2(0);
if (same_side(A, B, p2, s1)) {
if (!same_side(A, B, p1, s2))
return false;
if (s1 != s2 || (s1 * s2 == 0)) {
line = mp(A, B);
return true;
}
}
return false;
}
bool try_candidates(int index_1, int index_2, const polygon & p1, const polygon & p2, pair<point,point> & line){
int n_p1 = size(p1), n_p2 = size(p2);
fo (d1, -1, 1) fo (d2, -1, -1) {
int index_1a = (3 * n_p1 + d1 + index_1) % n_p1;
int index_1b = (index_1 + n_p1 + 1) % n_p1;
int index_2a = (3 * n_p2 + d2 + index_2) % n_p2;
int index_2b = (index_2 + n_p2 + 1) % n_p2;
if (are_equal(p1[index_1a], p2[index_2a]))
continue;
int s1 = 0, s2 = 0;
point A = p1[index_1a], B = p2[index_2a];
if (same_side(A, B, p2, s1)) {
if (!same_side(A, B, p1, s2))
continue;
if (s1 != s2 || (s1 * s2 ==0)) {
line = mp(A, B);
return true;
}
}
}
return false;
}
inline bool try_closest(int index1, int index2, const polygon & p1, const polygon & p2, pair<point, point> & line) {
double res1 = 1e20, res2 = 1e20;
int pair_1 = -1, pair_2 = -1;
int n_p1 = size(p1), n_p2 = size(p2);
re (i, size(p2)) {
double cur_dist = distance_point_segment(p1[index1], p2[i], p2[(i + 1) % size(p2)]);
if (cur_dist < res1) {
res1 = cur_dist;
pair_1 = i;
}
}
re (i, size(p1)) {
double cur_dist = distance_point_segment(p2[index2], p1[i], p1[(i+1) % size(p1)]);
if (cur_dist < res2) {
res2 = cur_dist;
pair_2 = i;
}
}
fo (d, -1, 1) {
if (try_pair(p1[(index1 + n_p1 + d) % n_p1], p1[(index1 + n_p1 + d + 1) % n_p1], p1, p2, line))
return true;
if (try_pair(p2[(index2 + n_p2 + d) % n_p2], p2[(index2 + n_p2 + d + 1) % n_p2], p1, p2, line))
return true;
if (try_pair(p1[(pair_2 + n_p1 + d) % n_p1], p1[(pair_2 + n_p1 + d + 1) % n_p1], p1, p2, line))
return true;
if (try_pair(p2[(pair_1 + n_p2 + d) % n_p2], p2[(pair_1 + n_p2 + d + 1) % n_p2], p1, p2, line))
return true;
}
/*
if (try_candidates(index1, pair_1, p1, p2, line))
return true;
if (try_candidates(pair_2, index2, p1, p2, line))
return true;
*/
return false;
}
/*
bool brute_test(const polygon & p1, const polygon & p2, pair<point,point> & line) {
bool ok = false;
re (i, size(p1))
if (try_candidates(i, p1, p2, line)) {
debug("on lower", p1[i], p1[(i + 1) % (size(p1))]);
debug("dists:", distance_point_polygon(p1[i], p2), distance_point_polygon(p1[(i + 1) % (size(p1))], p2));
//return true;
ok = true;
}
re (i, size(p2))
if (try_candidates(i, p2, p1, line)) {
debug("on upper", p2[i], p2[(i + 1) % (size(p2))]);
debug("dists:", distance_point_polygon(p2[i], p1), distance_point_polygon(p2[(i + 1) % (size(p2))], p1));
ok = true;
// return true;
}
assert(ok);
return false;
}
*/
pair<point, point> solve() {
if (size(lower) == 1) {
return mp(lower.front(), upper.front());
}
pair<point, point> line;
lower_hull = convex_hull(lower);
upper_hull = convex_hull(upper);
return line;
//debug("upper hull:", upper_hull);
//debug("lower hull:", lower_hull);
//brute_test(lower_hull, upper_hull, line);
//return line;
// find closest point in lower_hull relative to upper_hull
int closest_lower = closest_point_to_polygon(lower_hull, upper_hull);
int closest_upper = closest_point_to_polygon(upper_hull, lower_hull);
int n_lower = size(lower_hull);
int n_upper = size(upper_hull);
if (try_closest(closest_lower, closest_upper, lower_hull, upper_hull, line))
return line;
assert(size(lower) == 1 || size(upper) == 1);
return mp(lower_hull[rand() % n_lower], upper_hull[rand() % n_upper]);
}
int main() {
//FILE * fin = stdin;
//FILE * fin = fopen("oypara_test.in", "r");
FILE * fin = fopen(filein.c_str(), "r");
FILE * fout = fopen(fileout.c_str(), "w");
Sc.init(fin);
N = Sc.nextInt();
lower.reserve(N + 1);
upper.reserve(N + 1);
int last_x = -1, all_vertical = 1;
pair<point, point> line;
fo (i, 1, N) {
int x = Sc.nextInt(), y1 = Sc.nextInt(), y2 = Sc.nextInt();
lower.pb(point(x, y1));
upper.pb(point(x, y2));
if (i == 1) {
last_x = x;
} else {
if (x != last_x)
all_vertical = 0;
}
}
if (all_vertical) {
line = mp(point(last_x, 0), point(last_x, 1));
} else {
line = solve();
}
//debug("line:", line);
fprintf(fout, "%d %d %d %d\n", (int)(real(line.first) + EPSILON),
(int)(imag(line.first) + EPSILON),
(int)(real(line.second) + EPSILON),
(int)(imag(line.second) + EPSILON));
return 0;
}
Catalin Tiseanu Zeus