#include <bits/stdc++.h>
#define sz(X) (int)((X).size())
#define entire(X) X.begin(),X.end()
using namespace std; typedef long long ll;
void OJize(){cin.tie(NULL);ios_base::sync_with_stdio(false);}
template <class T1, class T2>ostream&operator<<(ostream &os,pair<T1,T2>const&x){os<<'('<<x.first<<", "<<x.second<<')';return os;}
template <class Ch, class Tr, class Container>basic_ostream<Ch,Tr>&operator<<(basic_ostream<Ch,Tr>&os,Container const&x){for(auto&y:x)os<<y<<" ";return os;}
typedef long double T;
// negative -1, zero 0, positive 1
template <typename T> int sgn(T x){return (T(0)<x) - (x<T(0));}
/////////////////////////////////// 2D POINT BASICS
struct pnt{
T x, y;
pnt operator+(pnt p){return {x+p.x, y+p.y};}
pnt operator-(pnt p){return {x-p.x, y-p.y};}
pnt operator*(T d){return {x*d, y*d};}
// !!!DOUBLE ONLY!!!
pnt operator/(T d){return {x/d, y/d};}
};
bool operator==(pnt a, pnt b){return a.x==b.x && a.y==b.y;}
bool operator!=(pnt a, pnt b){return !(a==b);}
bool operator<(pnt a, pnt b){return a.x<b.x || (a.x==b.x && a.y<b.y);}
bool operator>(pnt a, pnt b){return b<a;}
bool operator<=(pnt a, pnt b){return a.x<b.x || (a.x==b.x && a.y<=b.y);}
bool operator>=(pnt a, pnt b){return b<=a;}
T sq(pnt p){return p.x*p.x + p.y*p.y;}
double abs(pnt p){return sqrt(sq(p));}
T distsq(pnt p, pnt q){return sq(p-q);}
double dist(pnt p, pnt q){return abs(p-q);}
// dot/cross
T dot(pnt v, pnt w){return v.x*w.x + v.y*w.y;}
bool is_perp(pnt v, pnt w){return dot(v,w) == 0;}
double angle(pnt v, pnt w){
return acos(clamp(dot(v,w) / abs(v) / abs(w), -1., 1.));
}
T cross(pnt v, pnt w){return v.x*w.y - v.y*w.x;}
// ccw (negative cw, zero collinear, positive ccw)
T orient(pnt a, pnt b, pnt c){return cross(b-a, c-a);}
int ccw(pnt a, pnt b, pnt c){return sgn(orient(a,b,c));}
const pnt O = {0,0};
// IO
istream& operator>>(istream& is, pnt &p){
return is >> p.x >> p.y;
}
ostream& operator<<(ostream& os, pnt p){
return os << "("<<p.x<<", "<<p.y<<")";
}
/////////////////////////////////// 2D POINT TRANSFORM
// translate p by vector v
pnt translate(pnt p, pnt v){return p+v;}
// scale p, centered at c, by factor
pnt scale(pnt p, pnt c, T factor){return c + (p-c)*factor;}
// rotate p, centered at c, by rad ccw
// !!!DOUBLE ONLY!!!
pnt rotate(pnt p, pnt c, double rad){
pnt q = p-c;
return pnt{q.x*cos(rad)-q.y*sin(rad), q.x*sin(rad)+q.y*cos(rad)} + c;
}
// rotate p, centered at c, by pi/2 ccw
pnt rot90(pnt p, pnt c){
pnt q = p-c;
return pnt{-q.y, q.x} + c;
}
pnt rot90(pnt p){return {-p.y, p.x};}
// transform r with the linear transformation that moves p--q to fp--fq
// !!!DOUBLE ONLY!!!
pnt linear_trans(pnt p, pnt q, pnt r, pnt fp, pnt fq){
pnt pq = q-p, num{cross(pq, fq-fp), dot(pq, fq-fp)};
return fp + pnt{cross(r-p, num), dot(r-p, num)} / sq(pq);
}
/////////////////////////////////// 2D POINT SORTING
// sort points by (x, y)
void cartesian_sort(vector<pnt> &V){
sort(entire(V), [](pnt v, pnt w){
return make_pair(v.x, v.y) < make_pair(w.x, w.y);
});
}
// sort points by angle, ccw, centered at c
// tiebreaker is distance
bool pnt_half(pnt p){
assert(p.x != 0 || p.y != 0);
return p.y > 0 || (p.y == 0 && p.x < 0);
}
void polar_sort(vector<pnt> &V, pnt c){
for(pnt &p: V) p = p-c;
sort(entire(V), [](pnt v, pnt w){
return make_tuple(pnt_half(v), T(0), sq(v)) <
make_tuple(pnt_half(w), cross(v,w), sq(w));
});
for(pnt &p: V) p = p+c;
}
// sort points by orthogonal projection on v
void line_sort(vector<pnt> &V, pnt v){
sort(entire(V), [&](pnt p, pnt q){
return dot(v,p) < dot(v,q);
});
}
/////////////////////////////////// straight line
// [1] direction vector v, y-offset c (NOT Y-INTERCEPT!)
// [2] equation ax+by = c
// [3] two points p, q
struct line{
pnt v; T c;
line(): v({0,0}), c(0) {}
line(pnt V, T C): v(V), c(C) {}
line(T A, T B, T C): v({B,-A}), c(C) {}
line(pnt P, pnt Q): v(Q-P), c(cross(v,P)) {}
// positive ccw, zero collinear, negaitve cw
T side(pnt p){return cross(v,p)-c;}
double sqdist(pnt p){return side(p)*side(p) / sq(v);}
double dist(pnt p){return abs(side(p)) / abs(v);}
bool cmp_proj(pnt p, pnt q){return dot(v,p) < dot(v,q);}
};
// line perpendicular to L, through p
line perp_line(line L, pnt p){return {p, p+rot90(L.v)};}
// translate L by vector v
line translate(line L, pnt v){return {L.v, L.c+cross(L.v,v)};}
// translate L orthogonally left by d
// !!!DOUBLE ONLY!!!
line shift_left(line L, double d){return {L.v, L.c+d*abs(L.v)};}
// line intersection, true iff unique
pair<bool, pnt> line_inter(line L1, line L2){
T d = cross(L1.v, L2.v);
if(d == 0) return {false, {0,0}};
return {true, (L2.v*L1.c - L1.v*L2.c) / d};
}
/////////////////////////////////// straight line, more advanced
// projection of p onto L
// !!!DOUBLE ONLY!!!
pnt projection(line L, pnt p){return p - rot90(L.v)*L.side(p)/sq(L.v);}
// reflection of p by L
// !!!DOUBLE ONLY!!!
pnt reflection(line L, pnt p){return p - rot90(L.v)*2*L.side(p)/sq(L.v);}
// angle bisection
// interior = true gives a bisector between L1 and L2
// false gives the other one
// !!!DOUBLE ONLY!!!
line angle_bisector(line L1, line L2, bool interior){
assert(cross(L1.v, L2.v) != 0);
double sign = interior? 1: -1;
return {L2.v/abs(L2.v) + L1.v/abs(L1.v) * sign,
L2.c/abs(L2.v) + L1.c/abs(L1.v) * sign};
}
/////////////////////////////////// line segment a--b
// whether p is in the circle with diameter a--b
bool in_disk(pnt a, pnt b, pnt p){
return dot(a-p, b-p) <= 0;
}
// whether p is in a--b
bool on_seg(pnt a, pnt b, pnt p){
return orient(a,b,p) == 0 && in_disk(a,b,p);
}
// whether a--b and c--d has a unique proper intersection
pair<bool, pnt> seg_inter_proper(pnt a, pnt b, pnt c, pnt d){
T oa = orient(c,d,a), ob = orient(c,d,b);
T oc = orient(a,b,c), od = orient(a,b,d);
if(sgn(oa)*sgn(ob) >= 0 || sgn(oc)*sgn(od) >= 0)
return {false, {0,0}};
return {true, (a*ob - b*oa) / (ob-oa)};
}
// intersection points (or their endpoints, if infinite)
// !!!DOUBLE ONLY!!!
vector<pnt> seg_inters(pnt a, pnt b, pnt c, pnt d) {
auto [res, out] = seg_inter_proper(a,b,c,d);
if(res) return {out};
vector<pnt> s;
if (on_seg(c,d,a)) s.push_back(a);
if (on_seg(c,d,b)) s.push_back(b);
if (on_seg(a,b,c) && c!=a && c!=b) s.push_back(c);
if (on_seg(a,b,d) && d!=a && d!=b) s.push_back(d);
return s;
}
/////////////////////////////////// line segment distance
// distance from p to a--b
// !!!DOUBLE ONLY!!!
double seg_pnt_dist(pnt a, pnt b, pnt p){
if(a != b){
line L(a,b);
if(L.cmp_proj(a,p) && L.cmp_proj(p,b))
return L.dist(p);
}
return min(abs(p-a), abs(p-b));
}
// distance from a--b to c--d
// !!!DOUBLE ONLY!!!
double seg_seg_dist(pnt a, pnt b, pnt c, pnt d){
if(seg_inter_proper(a,b,c,d).first) return 0;
return min({seg_pnt_dist(a,b,c), seg_pnt_dist(a,b,d),
seg_pnt_dist(c,d,a), seg_pnt_dist(c,d,b)});
}
/////////////////////////////////// polygon
typedef vector<pnt> Polygon;
// twice the area
// positive ccw, negative cw
T area2(pnt a, pnt b, pnt c){return cross(b-a, c-a);}
T area2(Polygon &P){
T a = 0; int n = sz(P);
for(int i=0; i<n; i++){
a+= cross(P[i], P[(i+1)%n]);
}
return a;
}
// number of integer points on the line
// !!!INTEGER ONLY!!!
ll points_on_line(pnt p, pnt q){
ll dx = abs(p.x - q.x), dy = abs(p.y - q.y);
return gcd(dx, dy) + 1;
}
// number of boundary & interior points on the line
// area = i + p/2 - 1
// !!!INTEGER ONLY!!!
pair<ll, ll> picks_theorem(Polygon &P){
if(sz(P) == 1) return {1ll, 0ll};
ll perip = points_on_line(P[0], P.back()) - sz(P);
for(int i=1; i<sz(P); i++) perip+= points_on_line(P[i], P[i-1]);
ll A2 = abs(area2(P));
ll inp = (A2 + 2 - perip) / 2;
return {perip, inp};
}
// TODO point in polygon testing
/////////////////////////////////// convex polygon
struct ConvexPolygon{
Polygon pnts, lower, upper;
ConvexPolygon() {} // empty polygon
ConvexPolygon(Polygon P){
auto it = min_element(entire(P));
rotate(P.begin(), it, P.end());
auto uppercut = max_element(entire(P));
pnts = P;
lower.insert(lower.end(), P.begin(), uppercut+1);
upper.insert(upper.end(), uppercut, P.end());
upper.push_back(P[0]);
reverse(entire(upper));
}
bool empty(){return pnts.empty();}
};
T area2(ConvexPolygon &P){return area2(P.pnts);}
// io
ostream& operator<<(ostream& os, ConvexPolygon P){
for(pnt p: P.pnts) os << p << ' ';
return os;
}
ConvexPolygon convex_hull(Polygon P){
if(sz(P) == 1) return ConvexPolygon(P);
Polygon U, L;
sort(entire(P));
for(pnt q: P){
while(sz(U)>1 && ccw(*(U.end()-2), U.back(), q) >= 0) U.pop_back();
while(sz(L)>1 && ccw(*(L.end()-2), L.back(), q) <= 0) L.pop_back();
U.push_back(q), L.push_back(q);
}
reverse(entire(U));
L.insert(L.end(), U.begin()+1, U.end()-1);
return ConvexPolygon(L);
}
// for each point determine whether it's in P
// ASSUME PNT IS SORTED!! otherwise, set sorted = false
vector<bool> bulk_pnt_in_convex(ConvexPolygon P, vector<pnt> &ps, bool sorted=true){
if(!sorted) sort(entire(ps));
int il = 1, iu = 1;
vector<bool> res(sz(ps));
for(int i=0; i<sz(ps); i++){
pnt p = ps[i];
if(p < P.lower[0] || p > P.lower.back()){res[i] = false; continue;}
if(p == P.lower[0]){res[i] = true; continue;}
while(il < sz(P.lower)-1 && P.lower[il] <= p) il++;
while(iu < sz(P.upper)-1 && P.upper[iu] <= p) iu++;
bool ansl = orient(P.lower[il-1], p, P.lower[il]) <= 0;
bool ansu = orient(P.upper[iu-1], p, P.upper[iu]) >= 0;
res[i] = ansl & ansu;
}
return res;
}
// TODO point in convex testing in O(logn)
// pairs of {upper[i], lower[j]}
vector<pair<int, int>> rotating_calipers(ConvexPolygon &P){
vector<pair<int, int>> res;
vector<pnt> &U = P.upper, &L = P.lower;
int i = 0, j = sz(L)-1;
while(i < sz(U)-1 || j > 0){
res.push_back({i, j});
if(i == sz(U)-1) j--;
else if(j == 0) i++;
else if((U[i+1].y-U[i].y)*(L[j].x-L[j-1].x) >
(L[j].y-L[j-1].y)*(U[i+1].x-U[i].x)) i++;
else j--;
}
return res;
}
pair<pnt, pnt> farthest_convex(ConvexPolygon &P){
assert(!P.empty());
T ansd = 0;
pair<pnt, pnt> ans;
for(auto [i, j]: rotating_calipers(P)){
T d = distsq(P.upper[i], P.lower[j]);
if(ansd < d) ansd = d, ans = {P.upper[i], P.lower[j]};
}
return ans;
}
// each halfplane is the left (.side > 0) region of each line
// !!! ASSUME BOUNDED !!!
// !!! ASSUME UNIQUE SLOPES !!! todo remove this assumption
ConvexPolygon halfplane_intersection(vector<line> planes){
// TODO sort(P, ...);
deque<line> D;
auto bad = [&](line A, line B, line L){
auto [b, p] = line_inter(A, B);
return b && L.side(p) <= 0;
};
for(line &L: planes){
while(sz(D) >= 2 && bad(D[sz(D)-2], D.back(), L)) D.pop_back();
while(sz(D) >= 2 && bad(D[0], D[1], L)) D.pop_front();
if(sz(D) < 2 || !bad(D.back(), L, D[0])) D.push_back(L);
}
if(sz(D) < 3) return ConvexPolygon();
Polygon P;
for(int i=0; i<sz(D); i++){
auto [b, p] = line_inter(D[i], D[(i+1)%sz(D)]);
P.push_back(p);
}
return ConvexPolygon(P);
}
/////////////////////////////////// circle
struct circle{
pnt o; T r;
circle(pnt P, T R): o(P), r(R) {}
// -1 interior, 0 border, 1 exterior
int side(pnt p){return sgn(distsq(o,p)-r*r);}
};
bool operator==(circle a, circle b){return a.o==b.o && a.r==b.r;}
bool operator!=(circle a, circle b){return !(a==b);}
// !!!DOUBLE ONLY!!!
circle circumcircle(pnt a, pnt b, pnt c){
b = b-a, c = c-a;
assert(cross(b,c) != 0);
pnt cen = a + rot90(b*sq(c) - c*sq(b))/cross(b,c)/2;
return {cen, dist(cen, a)};
}
// circle-line intersection
// this only works for coords <= 30,000. TODO make it larger
// WARNING: UNTESTED!!!
int circle_line_inter0(circle C, line L){
pnt o = C.o; T r = C.r;
return 1 + sgn(r*r*sq(L.v) - L.side(o)*L.side(o));
}
// !!!DOUBLE ONLY!!!
// WARNING: UNTESTED!!!
tuple<int,pnt,pnt> circle_line_inter(circle C, line L){
pnt o = C.o; T r = C.r;
double h2 = r*r - L.sqdist(o);
if(h2 >= 0){
pnt p = projection(L, o), h = L.v*sqrt(h2) / abs(L.v);
return {1+sgn(h2), p-h, p+h};
}
return {0,{0,0},{0,0}};
}
// circle-seg intersection
// this only works for coords <= 30,000. TODO make it larger
// WARNING: UNTESTED!!!
int circle_seg_inter0(circle C, pnt a, pnt b){
T r = C.r; a = a-C.o, b = b-C.o;
pnt v = b-a;
T f0 = sq(a)-r*r, f1 = sq(b)-r*r, d = dot(a,v);
if(-sq(v) <= d && d <= 0){
T fm = sq(v)*f0 - d*d;
if(fm >= 0) return 1 - sgn(fm);
return (f0 >= 0) + (f1 >= 0);
}
return !((f0>0 && f1>0) || (f0<0 && f1<0));
}
// TODO make intersection finding version
// make use of circle line intersection and in_disk
// (not on_seg because of float errors!)
// circle-circle intersection
// -1 if infinite
// this only works for coords <= 30,000. TODO make it larger
int circle_circle_inter0(circle A, circle B){
if(A == B) return -1;
pnt o1 = A.o, o2 = B.o, d = B.o-A.o;
T r1 = A.r, r2 = B.r, d2 = sq(d);
if(d2 == 0){return 0;}
T R = d2+r1*r1-r2*r2;
return 1+sgn(4*r1*r1*d2 - R*R);
}
// !!!DOUBLE ONLY!!!
tuple<int,pnt,pnt> circle_circle_inter(circle A, circle B){
if(A == B) return {-1,{0,0},{0,0}};
pnt o1 = A.o, o2 = B.o, d = B.o-A.o;
T r1 = A.r, r2 = B.r, d2 = sq(d);
if(d2 == 0){return {0,{0,0},{0,0}};}
double pd = (d2 + r1*r1 - r2*r2)/2;
double h2 = r1*r1 - pd*pd/d2;
if(h2 >= 0){
pnt p = o1 + d*pd/d2, h = rot90(d)*sqrt(h2/d2);
return {1+sgn(h2), p-h, p+h};
}
return {0,{0,0},{0,0}};
}
// TODO tangent
///////////////////////////////////////////////
/////////// MAIN CODE BELOW ///////////////////
///////////////////////////////////////////////
int main(){OJize();
cout << fixed << setprecision(3);
int n; cin>>n;
Polygon P(n);
for(int i=0; i<n; i++) cin>>P[i];
vector<line> L(n);
for(int i=0; i<n; i++) L[i] = line(P[i], P[(i+1)%sz(P)]);
double lo = 0, hi = 40000000;
for(int ITER=0; ITER<60; ITER++){
double d = (lo + hi) / 2;
vector<line> LS(n);
for(int i=0; i<n; i++) LS[i] = shift_left(L[i], d);
ConvexPolygon C = halfplane_intersection(LS);
if(C.empty()){hi = d; continue;}
auto [a, b] = farthest_convex(C);
if(dist(a, b) < 2*d) hi = d;
else lo = d;
}
cout << lo;
}
Compilation message
twocircles.cpp: In function 'double angle(pnt, pnt)':
twocircles.cpp:40:58: error: no matching function for call to 'clamp(T, double, double)'
40 | return acos(clamp(dot(v,w) / abs(v) / abs(w), -1., 1.));
| ^
In file included from /usr/include/c++/10/algorithm:62,
from /usr/include/x86_64-linux-gnu/c++/10/bits/stdc++.h:65,
from twocircles.cpp:1:
/usr/include/c++/10/bits/stl_algo.h:3675:5: note: candidate: 'template<class _Tp> constexpr const _Tp& std::clamp(const _Tp&, const _Tp&, const _Tp&)'
3675 | clamp(const _Tp& __val, const _Tp& __lo, const _Tp& __hi)
| ^~~~~
/usr/include/c++/10/bits/stl_algo.h:3675:5: note: template argument deduction/substitution failed:
twocircles.cpp:40:58: note: deduced conflicting types for parameter 'const _Tp' ('long double' and 'double')
40 | return acos(clamp(dot(v,w) / abs(v) / abs(w), -1., 1.));
| ^
In file included from /usr/include/c++/10/algorithm:62,
from /usr/include/x86_64-linux-gnu/c++/10/bits/stdc++.h:65,
from twocircles.cpp:1:
/usr/include/c++/10/bits/stl_algo.h:3693:5: note: candidate: 'template<class _Tp, class _Compare> constexpr const _Tp& std::clamp(const _Tp&, const _Tp&, const _Tp&, _Compare)'
3693 | clamp(const _Tp& __val, const _Tp& __lo, const _Tp& __hi, _Compare __comp)
| ^~~~~
/usr/include/c++/10/bits/stl_algo.h:3693:5: note: template argument deduction/substitution failed:
twocircles.cpp:40:58: note: deduced conflicting types for parameter 'const _Tp' ('long double' and 'double')
40 | return acos(clamp(dot(v,w) / abs(v) / abs(w), -1., 1.));
| ^
