#1022173 (C++23(gcc13)) No.724 円と円の間の円

提出ソース
結果

問題	No.724 円と円の間の円
ユーザー	penguinman
提出日時	2024-10-25 01:17:02
言語	C++23(gcc13) (gcc 13.2.0 + boost 1.83.0)
結果	AC
実行時間	3 ms / 2,000 ms
コード長	15,136 bytes
コンパイル時間	5,701 ms
コンパイル使用メモリ	315,900 KB
実行使用メモリ	6,820 KB
最終ジャッジ日時	2024-10-25 01:17:09
合計ジャッジ時間	6,985 ms
ジャッジサーバーID （参考情報）	judge2 / judge3
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 6,816 KB
testcase_01	AC	2 ms 6,820 KB
testcase_02	AC	2 ms 6,820 KB
testcase_03	AC	3 ms 6,820 KB
testcase_04	AC	2 ms 6,816 KB
testcase_05	AC	2 ms 6,816 KB
testcase_06	AC	2 ms 6,820 KB
testcase_07	AC	2 ms 6,820 KB
testcase_08	AC	2 ms 6,820 KB
testcase_09	AC	2 ms 6,820 KB
testcase_10	AC	2 ms 6,816 KB
testcase_11	AC	2 ms 6,820 KB
testcase_12	AC	2 ms 6,816 KB
testcase_13	AC	2 ms 6,816 KB
testcase_14	AC	3 ms 6,820 KB
testcase_15	AC	2 ms 6,820 KB
testcase_16	AC	2 ms 6,820 KB
testcase_17	AC	2 ms 6,820 KB
testcase_18	AC	2 ms 6,820 KB
testcase_19	AC	2 ms 6,816 KB
testcase_20	AC	2 ms 6,820 KB
testcase_21	AC	2 ms 6,820 KB
testcase_22	AC	2 ms 6,820 KB
testcase_23	AC	2 ms 6,816 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <bits/stdc++.h>

#define rep(i, n) for(ll i = 0; i < ll(n); i++)
#define rep2(i, s, n) for(ll i = ll(s); i < ll(n); i++)
#define rrep(i, n) for(ll i = ll(n) - 1; i >= 0; i--)
#define rrep2(i, n, t) for(ll i = ll(n) - 1; i >= ll(t); i--)
#define all(a) a.begin(), a.end()
#define SZ(a) ll(a.size())
#define eb emplace_back
using namespace std;
using ll = long long;
using P = pair<ll, ll>;
using vl = vector<ll>;
using vvl = vector<vl>;
using vp = vector<P>;
using vvp = vector<vp>;
bool chmin(auto& a, auto b) { return a > b ? a = b, 1 : 0; }
bool chmax(auto& a, auto b) { return a < b ? a = b, 1 : 0; }

using ld = long double;
using Point = complex<ld>;
using std::norm;
using vd = vector<ld>;
using vdd = vector<vd>;
using pdd = pair<ld, ld>;

constexpr ld INF = 1e18;
constexpr ld EPS = 1e-6;

constexpr Point NULL_POINT = (INF, INF);
constexpr Point ORIGIN = (0, 0);

//reversed priority_queue
template<class T>
class prique :public std::priority_queue<T, std::vector<T>, std::greater<T>> {};

inline bool isNullPoint(Point &p){
    return p.real() > INF/2;
}

inline bool equal(const ld &a, const ld &b){
    return fabsl(a-b) < EPS;
}

inline bool equal(const Point &a, const Point &b){
    return equal(a.real(), b.real()) && equal(a.imag(), b.imag());
}


Point unitVector(const Point &a){
    return a / abs(a);
}

Point normalVector(const Point &a){
    return a * Point(0,1);
}

ld dot(const Point &a, const Point &b){
    return a.real() * b.real() + a.imag() * b.imag();
}

ld cross(const Point &a, const Point &b){
    return (a.real() * b.imag() - a.imag()*b.real());
}

Point rotate(const Point &p, const ld &theta){
    return p * Point(cos(theta), sin(theta));
}

ld radianToDegree(const ld &radian) { return radian * 180.0 / M_PI; }

ld degreeToRadian(const ld &degree) { return degree * M_PI / 180.0; }

struct Line {
    Point a, b;
    Line() = default;
    Line(Point a, Point b) : a(a), b(b) {}

    // Ax + By = C
    Line(ld A, ld B, ld C){
        if(equal(A, 0)){
            assert(!equal(B, 0));
            a = Point(0, C/B), b = Point(1, C/B);
        }
        else if(equal(B, 0)){
            a = Point(C/A, 0), b = Point(C/A, 1);
        }
        else{
            a = Point(C/A, 0), b = Point(0, C/B);
        }
    }
};

struct Segment : Line{
    Segment() = default;
    Segment(Point a, Point b) : Line(a, b){}
};

struct Circle {
    Point p;
    ld r;

    Circle() = default;

    Circle(Point p, ld r) : p(p), r(r) {}

    bool contain(const Point &q){
        return (abs(q-p) <= r+EPS);
    }
};

Point projection(const Line &l, const Point &p){
    ld t = dot(l.b-l.a, l.b-p) / norm(l.b-l.a);
    return l.a + (l.b-l.a) * t;
}

Point reflection(const Line &l, const Point &p){
    return p + (projection(l, p) - p) * ld(2.0);
}

// a → b → c の位置関係

/**
 * a → b → c の位置関係
 * 
 * @return
 * 1 : 反時計回り
 * -1 : 時計回り
 * 2 : c → a → b
 * -2 : a → b → c
 * 0 : a → c → b
 */
int ccw(const Point &a, Point b, Point c){
    b -= a, c -= a;
    if(cross(b,c) > EPS){
        return 1;
    }
    else if(cross(b,c) < -EPS){
        return -1;
    }
    else if(dot(b,c) < 0){
        return 2;
    }
    else if(norm(b) < norm(c)){
        return -2;
    }
    return 0;
}

// 直行？
bool isOrthogonal(const Line &a, const Line &b){
    return equal(dot(a.b-a.a, b.b-b.a), 0);
}

// 並行？
bool isParallel(const Line &a, const Line &b){
    return equal(cross(a.b-a.a, b.b-b.a), 0);
}

// 交差判定 → 線分 ab に対して c, d が反対側にあるか？
bool isIntersect(const Segment &s, const Segment &t){
    return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0
    && ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
}

//交点
Point crossPoint(const Line &s, const Line &t){
    ld d1 = cross(s.b-s.a, t.b-t.a);
    ld d2 = cross(s.b-s.a, s.b-t.a);
    if(equal(fabs(d1), 0)){
        if(equal(fabs(d2), 0)) return t.a;
        else return NULL_POINT;
    }
    return t.a + (t.b-t.a) * d2 / d1;
}

Point crossPoint(const Segment &s, const Segment &t){
    if(!isIntersect(s, t)) return NULL_POINT;
    return crossPoint(Line(s), Line(t));
}

ld distanceBetweenLineAndPoint(const Line &l, const Point &p){
    return fabs(cross(l.b-l.a, p-l.a)) / abs(l.b-l.a);
}

ld distanceBetweenSegmentAndPoint(const Segment &l, const Point &p){
    if(dot(l.b-l.a, p-l.a) < EPS) {
        return abs(p-l.a);
    }
    if(dot(l.a-l.b, p-l.b) < EPS){
        return abs(p-l.b);
    }
    return fabs(cross(l.b-l.a, p-l.a)) / abs(l.b-l.a);
}

ld distanceBetweenSegments(const Segment &s, const Segment &t){
    if(isIntersect(s, t)) return (ld)0;
    return std::min({
        distanceBetweenSegmentAndPoint(s, t.a),
        distanceBetweenSegmentAndPoint(s, t.b),
        distanceBetweenSegmentAndPoint(t, s.a),
        distanceBetweenSegmentAndPoint(t, s.b),
    });
}

// v : 多角形の頂点、反時計回り
struct Polygon{
    vector<Point> v;
    Polygon() = default;
    Polygon(vector<Point> p) : v(p) {}
    
    void add(Point a) {
        v.eb(a);
    }

    void rem(){
        assert(!v.empty());
        v.pop_back();
    }

    ld getArea(){
        ld res = 0;
        int n = v.size();
        assert(n >= 1);
        rep(i,n){
            res += cross(v[i], v[(i+1)%n]);
        }
        return res * ld(0.5);
    }

    /**
     * TODO : 一直線上にある場合の判定について少し考える
     */
    bool isConvex(){
        int n = v.size();
        rep(i,n){
            Point pre = v[i];
            Point now = v[(i+1)%n];
            Point nxt = v[(i+2)%n];
            if(ccw(pre, now, nxt) == -1) return false;
        }
        return true;
    }
    
    /**
     * @return
     * 2 : 含まれる
     * 1 : 辺上
     * 0 : 含まれない
     */
    int contain(Point &p){
        bool in = false;
        int n = v.size();
        rep(i,n){
            Point a = v[i]-p, b = v[(i+1)%n]-p;
            if(imag(a) > imag(b)) {
                std::swap(a, b);
            }
            if(imag(a) <= EPS && EPS < imag(b) && cross(a, b) < -EPS){
                in = !in;
            }
            if(equal(cross(a, b), 0) && dot(a, b) <= EPS){
                return 1;
            }
        }
        return (in ? 2 : 0);
    }
};

// 凸包、反時計回り
// @param border
// true -> border is included
// false -> border is not included
Polygon ConvexHull(Polygon p, bool border){
    ld feps = EPS;
    if(border) feps = -EPS;

    int n = p.v.size();

    if (n <= 1) return p;

    sort(all(p.v), [](const Point &a, const Point &b){
        if(equal(real(a), real(b))) return imag(a) < imag(b);
        return real(a) < real(b);
    });

    vector<Point> ret(2 * n);
    int k = 0;

    // 一直線上の3点を含める -> (< -EPS)
    // 含め無い -> (< EPS)
    rep(i,n){
        while(k >= 2){
            Point pre = ret[k-2], now = ret[k-1], nxt = p.v[i];
            if(cross(now-pre, nxt-now) < feps) k--;
            else break;
        }
        ret[k++] = p.v[i];
    }

    int t = k+1;

    rrep(i,n-1){
        while(k >= t){
            Point pre = ret[k-2], now = ret[k-1], nxt = p.v[i];
            if(cross(now-pre, nxt-now) < feps) k--;
            else break;
        }
        ret[k++] = p.v[i];
    }

    ret.resize(k-1);

    if (!border && ret.size() == 2 && equal(ret[0], ret[1])) ret.pop_back();
    
    return Polygon(ret);
}

/**
 * @return
 * 4 : 2 つの円が離れている
 * 3 : 外接している
 * 2 : 2 点で交わる
 * 1 : 内接している
 * 0 : 内包している
 */
int isIntersect(const Circle &c1, const Circle &c2){
    ld d = fabs(c1.p-c2.p);

    if(d > c1.r + c2.r + EPS) return 4;

    if(equal(d, c1.r + c2.r)) return 3;

    if(equal(d, fabs(c1.r-c2.r))) return 1;

    if(d < fabs(c1.r-c2.r)-EPS) return 0;

    return 2;
}

Circle inCircle(const Point &a, const Point &b, const Point &c){
    ld A = abs(a-c), B = abs(a-c), C = abs(a-b);

    Point p = a * A + b * B + c * C;
    p /= (A+B+C);

    ld r = distanceBetweenLineAndPoint(Line(a,b), p);

    return Circle(p, r);
}

vector<Point> crossPoint(const Circle &c, const Line &l){
    vector<Point> res;
    ld d = distanceBetweenLineAndPoint(l, c.p);

    if(d > c.r + EPS){
        return res;
    }
    
    Point h = projection(l, c.p);
    if(equal(d, c.r)){
        res.eb(h);
        return res;
    }

    Point e = unitVector(l.b-l.a);

    ld ph = sqrtl(c.r*c.r-d*d);

    res.eb(h-e*ph);
    res.eb(h+e*ph);

    return res;
}

vector<Point> crossPoint(const Circle &c1, const Circle &c2){
    vector<Point> res;
    int mode = isIntersect(c1, c2);

    ld d = abs(c1.p-c2.p);

    if(mode == 4 || mode == 0){
        return res;
    }

    if(mode == 3){
        ld t = c1.r/(c1.r+c2.r);
        res.eb(c1.p+(c2.p-c1.p)*t);
        return res;
    }

    if(mode == 1){
        if(c2.r < c1.r-EPS){
            res.eb(c1.p+(c2.p-c1.p)*(c1.r/d));
        }
        else{
            res.eb(c2.p+(c1.p-c2.p)*(c2.r/d));
        }
        return res;
    }

    // 2 円が重なる場合
    ld rc1 = (c1.r * c1.r + d * d - c2.r * c2.r) / (2 * d);
    ld rs1 = sqrtl(c1.r * c1.r - rc1 * rc1);

    if(c1.r - fabs(rc1) < EPS){
        rs1 = 0;
    }

    Point e12 = unitVector(c2.p-c1.p);

    res.eb(c1.p + rc1*e12 + rs1*normalVector(e12));
    res.eb(c1.p + rc1*e12 - rs1*normalVector(e12));

    return res;
}

vector<Point> tangentToCircle(const Point &p, const Circle &c){
    if(c.r < abs(c.p-p)-EPS) return {};
    if(equal(c.r, abs(c.p-p))) return {p};

    return crossPoint(c, Circle(p, sqrtl(norm(c.p-p)-c.r*c.r)));
}

/**
 * TODO : 完全理解
 */
vector<Line> tangent(const Circle &a, const Circle &b){
    vector<Line> ret;

    ld g = abs(a.p-b.p);

    if(equal(g, 0)){
        return ret;
    }

    Point u = unitVector(b.p-a.p);
    Point v = normalVector(u);

    for(int s : {-1, 1}){
        ld h = (a.r + b.r * s) / g;
        if(equal(h*h, 1)){
            ret.eb(a.p+(h>0 ? u : -u)*a.r, a.p+(h>0 ? u:-u)*a.r+v);
        }
        else if(1-h*h >= EPS){
            Point U = u*h, V = v*sqrtl(1-h*h);
            ret.eb(a.p+(U+V)*a.r, b.p-(U+V)*(b.r*s));
            ret.eb(a.p+(U-V)*a.r, b.p-(U-V)*(b.r*s));
        }
    }
    return ret;
}

ll mostCommonPoint(vector<Circle> C){
    ll N = C.size();
    vector<Point> cand;
    rep(i,N){
        cand.eb(C[i].p);
        rep(j,i){
            for(auto el: crossPoint(C[i], C[j])) cand.eb(el);
        }
    }
    ll ans = 0;
    for(auto center: cand){
        ll cnt = 0;
        rep(i,N){
            cnt += C[i].contain(center);
        }
        ans = std::max(ans, cnt);
    }
    return ans;
}

vector<Point> argSort(vector<Point> p) {
    sort(all(p), [](const Point &a, const Point &b){
        auto type = [&](Point x) {
            if (equal(x, ORIGIN)) return 0;
            if (x.imag() < -EPS || (equal(x.imag(), 0) && -EPS < x.real())) return -1;
            return 1;
        };
        int ta = type(a), tb = type(b);
        if (ta != tb) return (ta < tb);
        return cross(a, b) > 0;
    });
    return p;
}

Point input() {
    ld x, y; cin >> x >> y;
    return Point(x, y);
}

pair<Point, Point> closest_pair(vector<Point> &x, int left = 0, int right = -1){
    if (right == -1) {
        right = x.size();
        sort(all(x), [](const Point &a, const Point &b){
            if (equal(a.real(), b.real())) return a.imag() < b.imag();
            return a.real() < b.real();
        });
    }
    if (right-left <= 1) {
        return make_pair(Point(INF, INF), Point(-INF, -INF));
    }

    int mid = (left+right)/2;
    auto p = x[mid];
    pair<Point, Point> res;
    {
        auto l = closest_pair(x, left, mid);
        auto r = closest_pair(x, mid, right);
        if (abs(l.first-l.second) < abs(r.first-r.second)) res = l;
        else res = r;
    }

    auto within = [&](ld a, ld b) {
        return fabs(a-b) < abs(res.first-res.second)-EPS;
    };

    auto chminP = [&](Point a, Point b) {
        if (abs(a-b) < abs(res.first-res.second)) res = make_pair(a, b);
    };

    vector<Point> y, z;
    {
        int l = left, r = mid;
        while(l < mid && r < right) {
            if (x[l].imag() < x[r].imag()) {
                y.eb(x[l++]);
            } else {
                y.eb(x[r++]);
            }
        }
        while (l < mid) y.eb(x[l++]);
        while (r < right) y.eb(x[r++]);

        rep2(i,left,right) {
            x[i] = y[i-left];
            if (within(x[i].real(), p.real())){
                z.eb(x[i]);
            }
        }
    }

    rep(i,z.size()) {
        rep2(j,i+1,z.size()) {
            if (within(z[i].imag(), z[j].imag())) chminP(z[i], z[j]);
            else break;
        }
        rrep2(j,i,0) {
            if (within(z[i].imag(), z[j].imag())) chminP(z[i], z[j]);
            else break;
        }
    }
    return res;
}

pair<Point, Point> furthest_pair(vector<Point> x) {

    vector<Point> p = ConvexHull(Polygon(x), false).v;

    if(p.size() == 1) return make_pair(p[0], p[0]);
    if(p.size() == 2) return make_pair(p[0], p[1]);

    auto res = make_pair(p[0], p[1]);

    auto chmaxP = [&](Point a, Point b) {
        if (abs(a-b) > abs(res.first-res.second)) res = make_pair(a, b);
    };

    int now = 0;
    rep(i,p.size()) {
        while(true){
            chmaxP(p[i], p[now]);
            if(cross(p[(i+1)%p.size()]-p[i], p[(now+1)%p.size()]-p[now]) < -EPS) break;
            now++;
            now %= SZ(p);
        }
    }
    return res;
}

Point outerCenter(Point a, Point b, Point c){
    ld A = norm(b-c), B = norm(c-a), C = norm(a-b);
    ld S = cross(c-a, b-a);
    return (A*(B+C-A)*a + B*(C+A-B)*b + C*(A+B-C)*c) / (ld(4)*S*S);
}

// 最小包含円、期待値 O(N)
ld min_ball(vector<Point> p){
    random_device rnd;
    mt19937 mt(rnd());
    shuffle(all(p), mt);

    ll N = SZ(p);

    Circle c;

    auto upd = [&](int i, int j){
        c = Circle((p[i]+p[j])/ld(2), abs(p[i]-p[j])/ld(2));
    };

    auto upd2 = [&](int i, int j, int k){
        Point cent = outerCenter(p[i], p[j], p[k]);
        c = Circle(cent, abs(cent-p[i]));
    };

    if(N == 1) return ld(0);
    upd(0,1);

    rep2(i,2,N){
        if(c.contain(p[i])) continue;
        upd(0,i);
        rep(j,i){
            if(c.contain(p[j])) continue;
            upd(j,i);
            rep(k,j){
                if(c.contain(p[k])) continue;
                upd2(i,j,k);
            }
        }
    }
    return c.r;
}

int main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
    cout << std::fixed << std::setprecision(15);
    ll N;
    ld A,B; cin >> N >> A >> B;
    A *= 2;
    B *= 2;
    // sqrt(AB) による反転
    Point c;
    if(N%2){
        c = Point((A+B)/ld(2), ld(N/2)*(B-A));
    } else{
        c = Point((A+B)/ld(2), ld(N/2)*(B-A)-(B-A)/ld(2));
    }
    Point x = c, y = c, z = c;
    x = Point(A, imag(c));
    y = Point(B, imag(c));
    z = c - Point(0, (B-A)/ld(2));
    x *= A*B/norm(x);
    y *= A*B/norm(y);
    z *= A*B/norm(z);
    c = outerCenter(x,y,z);
    cout << abs(c-x) << endl;
}
yukicoder

結果

テストケース

ソースコード
