Submission #479338

#	Time	Username	Problem	Language	Result	Execution time	Memory
479338		jh05013	2circles (balkan11_2circles)	C++17	80 / 100	220 ms	7916 KiB

twocircles

#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 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 (stderr)

twocircles.cpp: In function 'int circle_circle_inter0(circle, circle)':
twocircles.cpp:432:9: warning: variable 'o1' set but not used [-Wunused-but-set-variable]
  432 |     pnt o1 = A.o, o2 = B.o, d = B.o-A.o;
      |         ^~
twocircles.cpp:432:19: warning: variable 'o2' set but not used [-Wunused-but-set-variable]
  432 |     pnt o1 = A.o, o2 = B.o, d = B.o-A.o;
      |                   ^~
twocircles.cpp: In function 'std::tuple<int, pnt, pnt> circle_circle_inter(circle, circle)':
twocircles.cpp:441:19: warning: variable 'o2' set but not used [-Wunused-but-set-variable]
  441 |     pnt o1 = A.o, o2 = B.o, d = B.o-A.o;
      |                   ^~

• Subtask 1: 80 / 100