Submission #638560

---

Source Code Expand

#include <bits/stdc++.h>//{{{
#define rep(i, n) for(int i = 0; i < (int)(n); ++i)
#define all(v) std::begin(v), std::end(v)
using namespace std;
template<class T>bool chmin(T&a,const T&b){if(a<=b)return false;a=b;return true;}
template<class T>bool chmax(T&a,const T&b){if(a>=b)return false;a=b;return true;}
//}}}

using Z = int;
using R = long double;
using P = complex<R>;
#define X real()
#define Y imag()
const R EPS = 1E-13;
const R PI = acos(R(-1));
int sgn(const R& r){ return (r > EPS) - (r < -EPS); }
int sgn(const R& a, const R &b){ return sgn(a - b); }
R dot(const P& a, const P& b){ return real(conj(a)*b); }
R det(const P& a, const P& b){ return imag(conj(a)*b); }
R angle(const P &p){ R theta = arg(p); return sgn(theta - PI) < 0 ? theta : -PI; }
R dPP(const P &a, const P &b){//{{{
    return abs(a - b);
}//}}}
struct L{//{{{
    P s, t;
    L(){}
    L(const P &s, const P &t) : s(s), t(t){}
};//}}}
struct C{//{{{
    P o; R r;
    C(){}
    C(const P &o, const R &r) : o(o), r(r){}
};//}}}
istream &operator>>(istream &is, P &p){//{{{
    R x, y; is >> x >> y;
    p.real(x); p.imag(y);
    return is;
}//}}}
namespace std{
    bool operator < (const P& a, const P& b) {//{{{
        return sgn(a.X, b.X) ? a.X < b.X : sgn(a.Y, b.Y) < 0;
    }//}}}
}
enum CCW{FRONT = 1, RIGHT = 2, BACK = 4, LEFT = 8, ON = 16, };
inline int ccw(const L& l, const P& p){//{{{
    P b = p - l.s, c = l.t - l.s;
    int scr = sgn(det(c, b));
    if(scr) return scr < 0 ? RIGHT : LEFT;
    if(sgn(dot(c, b)) < 0) return BACK;
    if(sgn(norm(c) - norm(b)) < 0) return FRONT;
    return ON;
}//}}}

struct Witch{//{{{
    P o;
    Z r, v;
    Witch(){}
    Witch(const P &o, const Z &r, const Z &v) : o(o), r(r), v(v){}
    friend istream &operator>>(istream &is, Witch &w){
        return is >> w.o >> w.r >> w.v;
    }
};//}}}

R dLP(const L& l, const P& p){//{{{
    return abs(det(l.t - l.s, p - l.s)) / abs(l.t - l.s);
}//}}}
// 二円の交点
vector<P> pS1S1(const C &a, const C &b){//{{{
    R d = dPP(a.o, b.o);
    R rc = (d * d + a.r * a.r - b.r * b.r) / (2 * d);
    if(a.r * a.r - rc * rc < 0) return {};
    R rs = sqrt(max(a.r * a.r - rc * rc, R(0)));
    P dir = (b.o - a.o) / d;
    return { a.o + dir * P(rc, rs), a.o + dir * P(rc, -rs) };
}//}}}
// 一点を通る, 円との接線
vector<P> tS1P(const C &c, const P &p){//{{{
    R d = dPP(c.o, p);
    if(d * d - c.r * c.r < 0) return {};
    R r = sqrt(max(d * d - c.r * c.r, R(0)));
    return pS1S1(c, C(p, r));
}//}}}
// 二円の共通接線の接点
vector<P> lS1S1(C a, C b){//{{{
    vector<P> res;
    if(a.r < b.r) swap(a, b);
    R d = dPP(a.o, b.o);
    P v = (b.o - a.o) / d;
    int s = sgn(a.r - b.r - d);
    if(s > 0) return res; // a が b を含む
    if(s == 0){ // a の内側に b が接している
        P t = a.o + a.r * v, u = t + v * P(0, 1);
        res.emplace_back(t);
        res.emplace_back(u);
        return res;
    }
    // 外接線
    if(sgn(a.r - b.r) == 0){
        P t = a.o + v * P(0, 1) * a.r;
        P u = a.o - v * P(0, 1) * a.r;
        res.emplace_back(t);
        res.emplace_back(t + b.o - a.o);
        res.emplace_back(u);
        res.emplace_back(u + b.o - a.o);
    }else{
        R theta = PI/2 - asin((a.r - b.r) / d);
        P e = polar(R(1), theta);
        res.emplace_back(a.o + v * a.r * e);
        res.emplace_back(b.o + v * b.r * e);
        res.emplace_back(a.o + v * a.r / e);
        res.emplace_back(b.o + v * b.r / e);
    }
    // 内接線
    s = sgn(d - a.r - b.r);
    if(s == 0){
        P t = a.o + a.r * v;
        res.emplace_back(t);
        res.emplace_back(t + v * P(0, 1));
    }else if(s > 0){
        R theta = PI/2 - asin((a.r + b.r) / d);
        P e = polar(R(1), theta);
        res.emplace_back(a.o + v * a.r * e);
        res.emplace_back(b.o - v * b.r * e);
        res.emplace_back(a.o + v * a.r / e);
        res.emplace_back(b.o - v * b.r / e);
    }
    return res;
}//}}}

R solve_easy(const vector<C> &cs, const P start, const P goal){//{{{
    auto is_valid_point = [&](const P &p){//{{{
        for(auto &c : cs) if(sgn(dPP(c.o, p), c.r) < 0) return false;
        return true;
    };//}}}
    auto is_valid_line = [&](const L &l){//{{{
        for(auto &c : cs){
            if(sgn(dLP(l, c.o), c.r) > 0) continue;
            const L m(c.o, c.o + (l.t - l.s) * P(0, 1));
            if(sgn(dPP(c.o, l.s), c.r) < 0) return false;
            if(sgn(dPP(c.o, l.t), c.r) < 0) return false;
            if((ccw(m, l.s) | ccw(m, l.t)) == (LEFT | RIGHT))
                return false;
        }
        return true;
    };//}}}
    auto is_valid_arc = [&](const C &c, const R a, const R b){//{{{
        P pa = c.o + polar(c.r, a), pb = c.o + polar(c.r, b);
        if(!is_valid_point(pa) or !is_valid_point(pb)) return false;
        P dir = (pa - c.o) + (pb - c.o);
        dir /= abs(dir);
        return is_valid_point(c.o + dir * c.r);
    };//}}}


    set<P> ps = {start, goal};
    for(auto &c : cs) for(auto &p : tS1P(c, start))
        ps.emplace(p);
    for(auto &c : cs) for(auto &p : tS1P(c, goal))
        ps.emplace(p);
    for(auto &a : cs) rep(dir, 4)
        ps.emplace(a.o + polar(a.r, PI / 2 * dir));
    for(auto &a : cs) for(auto &b : cs) for(auto &p : pS1S1(a, b))
        ps.emplace(p);
    for(auto &a : cs) for(auto &b : cs) for(auto &p : lS1S1(a, b))
        ps.emplace(p);
    // make graph
    map<P, map<P, R>> g;
    for(auto &p : ps) for(auto &q : ps) if(is_valid_line(L(p, q))){
        if(g[p].count(q)) chmin(g[p][q], dPP(p, q));
        else              g[p][q] = dPP(p, q);
    }
    for(auto &c : cs){//{{{
        vector<pair<R, P>> as;
        for(auto &p : ps) if(sgn(dPP(p, c.o), c.r) == 0)
            as.emplace_back(angle(p - c.o), p);
        sort(all(as));
        rep(i, as.size()){
            const int j = (i + 1) % as.size();
            R ai = as[i].first, aj = as[j].first;
            P pi = as[i].second, pj = as[j].second;
            if(is_valid_arc(c, ai, aj)){
                R a = aj - ai;
                if(a < 0) a += 2 * PI;
                if(g[pi].count(pj)) chmin(g[pi][pj], c.r * a);
                else                g[pi][pj] = c.r * a;
                if(g[pj].count(pi)) chmin(g[pj][pi], c.r * a);
                else                g[pj][pi] = c.r * a;
            }
        }
    }//}}}
    const R inf = 1E100;
    map<P, R> d;
    for(auto &p : ps) d[p] = inf;
    priority_queue<pair<R, P>, vector<pair<R, P>>, greater<pair<R, P>>> pq;
    pq.emplace(d[start] = 0, start);
    while(!pq.empty()){
        R c; P p; tie(c, p) = pq.top(); pq.pop();
        if(d[p] < c) continue;
        for(auto &qc : g[p]){
            P q; R cc; tie(q, cc) = qc;
            if(chmin(d[q], c + cc)) pq.emplace(d[q], q);
        }
    }
    return d[goal];
}//}}}

template<class R, class CMP = less<R>, class F>
R find_max(R l, R r, F f){//{{{
    constexpr R EPS = 1E-11;
    const R phi = R(2) / (3 + sqrt(R(5)));
    R ml = l + phi * (r - l), mr = r - phi * (r - l);
    R fml = f(ml), fmr = f(mr);
    while(mr - ml > EPS){
        if(fml < fmr){ // '<' : max, '>' : min
            l = ml; ml = mr; mr = r - phi * (r - l);
            fml = fmr; fmr = f(mr);
        }else{
            r = mr; mr = ml; ml = l + phi * (r - l);
            fmr = fml; fml = f(ml);
        }
    }
    return ml;
}//}}}

void solve_hard(Witch w, P start, P goal, R speed){//{{{
    start -= w.o;
    goal  -= w.o;
    w.o    = P(0, 0);
    auto rot = [&](P p){//{{{
        p /= abs(p);
        start /= p;
        goal  /= p;
    };//}}}
    rot(start);
    if(sgn(goal.Y) < 0) goal = conj(goal);
    // w.o は原点, start は実軸上, goal は上半平面.

    {   // 直接行く//{{{
        P dir = goal - start;
        //   goal が殺される                    goal に直接つく
        if(sgn((abs(goal) - w.r) / w.v, abs(dir) / speed) <= 0){
            cout << "Impossible" << endl;
            return;
        }
        auto f = [&](const R x){
            P p = x * start + (1-x) * goal;
            R t = abs(start - p) / speed;
            return (w.r + w.v * t) - abs(p);
        };
        if(sgn(f(find_max<R>(0, 1, f))) <= 0){
            cout << abs(dir) / speed << endl;
            return;
        }
    }//}}}
    auto dbl = [](const R &x){ return x * x; };

    const R t1 = sqrt((dbl(start.X) - dbl(w.r)) / (dbl(speed) - dbl(w.v)));
    if(sgn(w.r + t1 * w.v, abs(goal)) >= 0){
        cout << "Impossible" << endl;
        return;
    }

    P angle;
    {//{{{
        R xvdx = w.r * w.v + sqrt((dbl(speed) - dbl(w.v)) * (dbl(start.X) - dbl(w.r)));
        R xvdy = sqrt(dbl(start.X * speed) - xvdx);

        auto ps = pS1S1(C(w.o, w.r + t1 * w.v), C(start, t1 * speed));
        assert(!ps.empty());
        P p = ps.front().Y < 0 ? ps.back() : ps.front();
        // => p * amp ----> goal
        rot(p);

        p = P(abs(p), 0);
        angle = p - start;
    }//}}}
    // alpha * beta^theta
    const R alpha = w.r + t1 * w.v;
    const R beta = exp(w.v / sqrt(dbl(speed) - dbl(w.v)));
    auto f = [&](const R theta){//{{{
        P p = polar(alpha * pow(beta, theta), theta);
        P dir = p * angle;
        if(sgn(det(goal - p, p)) > 0) return R(+1);
        return det(goal - p, dir);
    };//}}}
    assert(f(0) < 0);
    assert(f(PI) > 0);
    R theta = [&]{//{{{
        R l = 0, r = PI;
        rep(_, 1000){
            R m = (l + r) / 2;
            (f(m) < 0 ? l : r) = m;
        }
        return l;
    }();//}}}
    const P p = polar(alpha * pow(beta, theta), theta);
    const R t12 = (abs(p) - w.r) / w.v;
    const R t3 = abs(p - goal) / speed;
    cout << t12 + t3 << endl;
}//}}}

int main(){
    cout << std::fixed << std::setprecision(10);

    P start, goal;
    R speed;
    int n;
    cin >> start >> goal >> speed >> n;
    // static witch, dymamic witch
    vector<C> witch_s;
    vector<Witch> witch_d;
    rep(_, n){
        Witch w; cin >> w;
        if(w.v == 0) witch_s.emplace_back(w.o, w.r);
        else         witch_d.emplace_back(w);
    }
    if(n == 1 and !witch_d.empty() and witch_d[0].v < speed){
        // hard problem
        solve_hard(witch_d.front(), start, goal, speed);
    }else{
        // easy problem
        const R res = solve_easy(witch_s, start, goal) / speed;
        if(res > 1E50){
            cout << "Impossible" << endl;
            return 0;
        }
        for(auto &w : witch_d) if(sgn((dPP(goal, w.o) - w.r) / w.v, res) < 0){
            cout << "Impossible" << endl;
            return 0;
        }
        cout << res << endl;
    }
    return 0;
}
// vim:set foldmethod=marker commentstring=//%s:

Submission Info

Judge Result