ASC 26 Problem I: Shepherd's Problem

#include <iostream>
#include <cstdio>
#include <iomanip>
#include <complex>
#include <vector>
#include <algorithm>
using namespace std;
typedef long double ld;
typedef complex<ld> pt;
const ld PI = acos(-1.L), EPS = 1e-12L;
ld cp(pt a, pt b) { return a.real()*b.imag()-b.real()*a.imag(); }
ld dp(pt a, pt b) { return a.real()*b.real()-a.imag()*b.imag(); }
// p = input points, a and b are global variables for results of intersections
vector<pt> p; pt a, b;
bool cc_inter(pt p1, int r1, pt p2, int r2) { // circle circle
    ld dq=norm(p1-p2), rq=r1*r1-r2*r2;
    ld dt=2*dq*(r1*r1+r2*r2)-dq*dq-rq*rq;
    a=b=(p1+p2)/2.L+(p2-p1)/2.L*rq/dq;
    if(dt<-EPS) return false;
    if(dt<EPS) return true;
    dt=sqrt(dt)/2/dq;
    a=b+(p1-p2)*pt(0,-1)*dt;
    b=b-(p1-p2)*pt(0,-1)*dt;
    return true;
}
bool cl_inter(pt c, int r, pt q1, pt q2) { // circle line
    pt u = q2-q1;
    ld uq=norm(u), d=uq*r*r-cp(u,q1-c)*cp(u,q1-c);
    if(d<-EPS) return false;
    ld s=dp(u, c-q1)/uq;
    a = b = q1+s*u;
    if(d<EPS) return true;
    d=sqrt(d)/uq;
    a = q1 + (s-d)*u, b = q1 + (s+d)*u;
    return true;
}
int main() {
    freopen("shepherd.in", "r", stdin);
    freopen("shepherd.out", "w", stdout);
    int n, d, D; cin >> n >> d >> D;
    for (int i=0; i < n; i++) {
        ld x, y; cin >> x >> y;
        p.push_back(pt(x,y));
    }

    // best = smallest angle, bestpt = point with smallest angle
    ld best=4; pt bestpt;
    // Test for every p[i] the candidate points on circle of radius D centered on p[i]
    for (int i=0; i < p.size(); i++) {
        // Endpoints of segment
        vector<ld> v; v.push_back(0); v.push_back(2*PI);
        // for every other point
        for (int j=0; j < p.size(); j++) if (j != i) {
            // small circle intersection
            if (cc_inter(p[i], D, p[j], d))
                v.push_back(arg(a-p[i])), v.push_back(arg(b-p[i]));
            // big circle intersection
            if (cc_inter(p[i], D, p[j], D))
                v.push_back(arg(a-p[i])), v.push_back(arg(b-p[i]));
            // line intersection with every pair of points
            for (int k=0; k < p.size(); k++) if (k != j) {
                if (cl_inter(p[i], D, p[j], p[k]))
                   v.push_back(arg(a-p[i])), v.push_back(arg(b-p[i]));
            }
        }
        // get rid of negative values and sort
        for (int j=0; j < v.size(); j++) if (v[j] < 0) v[j] += 2*PI;
        sort(v.begin(), v.end());

        // For each consecutive segment of candidate points
        for (int j=0, k=1; k < v.size(); j=k++) {
            ld lo = v[j], hi = v[k], m = (lo+hi)/2;
            // Either the whole segment is valid, or it's invalid, so just test the midpoint
            pt pp = p[i] + polar(D+0.L, m);
            bool valid = true;
            // For each point
            for (int r = 0; r < p.size(); r++) {
                // distance is valid?
                if (abs(pp-p[r]) + EPS < d || abs(pp-p[r]) - EPS > D) valid = false;
                // for each triangle, is the point contained in the triangle?
                for (int s = r+1; s < p.size(); s++) for (int t = s+1; t < p.size(); t++) {
                    ld dir = cp(p[s]-p[r], p[t]-p[r]);
                    if (cp(p[s]-p[r], pp-p[r]) * dir > EPS &&
                            cp(p[t]-p[s], pp-p[s]) * dir > EPS &&
                            cp(p[r]-p[t], pp-p[t]) * dir > EPS)
                        valid = false;
                }
            }
            if (!valid) continue;

            // The segment is valid, so ternary search on the segment
            while(true) {
                ld c1 = (lo*2+hi)/3, c2 = (lo+hi*2)/3;
                // p1 and p2 are candidate points to test
                pt p1 = p[i] + polar(D+0.L, c1), p2 = p[i] + polar(D+0.L, c2);
                ld a1=0, a2=0; // find angles for p1 and p2 (max over all pairs of points)
                for (int r=0; r < p.size(); r++) for (int s=r+1; s < p.size(); s++) {
                    ld t1 = abs(arg(p[r]-p1)-arg(p[s]-p1));
                    a1 = max(a1, min(t1, 2*PI-t1));
                    ld t2 = abs(arg(p[r]-p2)-arg(p[s]-p2));
                    a2 = max(a2, min(t2, 2*PI-t2));
                }
                if (a1 < a2) {
                    hi = c2;
                    if (a1 < best) best = a1, bestpt = p1;
                }
                else {
                    lo = c1;
                    if (a2 < best) best = a2, bestpt = p2;
                }
                if (hi-lo < EPS) break;
            }
        }
    }
    if (best > PI) cout << "impossible" << endl;
    else cout << fixed << setprecision(10) << bestpt.real() << " " << bestpt.imag() << endl;
    return 0;
}