2023-08-12

AtCoder ABC 314 Ex - Disk and Segments (赤色, 625 点)

AtCoder ABC-Ex 計算幾何円線分の交差判定や距離計算三分探索高度典型連続最適化連続量問題最小コスト凸関数山登り法二次元平面上のN点の問題 AtCoder625点赤色diff 最適化問題浮動小数点型を扱う問題

本当に三分探索で通るのね！　これらの問題を思い出した

drken1215.hatenablog.com

問題へのリンク

問題概要

二次元平面上に $N$ 本の線分がある。この平面上の円盤 (円周と内部を含む) のうち、 $N$ 本の線分すべてと共有点をもつようなものの、最小の半径を求めよ。

制約

$2 \le N \le 100$
線分の端点の x 座標や y 座標は $0$ 以上 $1000$ 以下

解法 (解説 AC)

点 $(x, y)$ を中心とした円盤のうち、条件を満たすものの半径の最小値を $f(x, y)$ としたとき、 $f(x, y)$ を最小化する問題と言える。

なお、 $f(x, y)$ の計算自体は「点 $(x, y)$ と $N$ 本の線分との距離の最大値」でよい。

$f(x, y)$ は凸関数になっているので、三分探索で解ける。

コード

#include <bits/stdc++.h>
using namespace std;

// Geometry basic settings
using DD = long double;
constexpr long double PI = 3.141592653589793238462643383279502884L;
constexpr long double INF = 1LL<<60;  // to be set appropriately
constexpr long double EPS = 1e-10;    // to be set appropriately
long double torad(int deg) {return (long double)(deg) * PI / 180;}
long double todeg(long double ang) {return ang * 180 / PI;}

// Point or Vector
struct Point {
    DD x, y;
    
    // constructor
    constexpr Point() : x(0), y(0) {}
    constexpr Point(DD x, DD y) : x(x), y(y) {}
    
    // various functions
    constexpr Point conj() const {return Point(x, -y);}
    constexpr DD dot(const Point &r) const {return x * r.x + y * r.y;}
    constexpr DD cross(const Point &r) const {return x * r.y - y * r.x;}
    constexpr DD norm() const {return dot(*this);}
    constexpr long double abs() const {return sqrt(norm());}
    constexpr long double amp() const {
        long double res = atan2(y, x);
        if (res < 0) res += PI*2;
        return res;
    }
    constexpr bool eq(const Point &r) const {return (*this - r).abs() <= EPS;}
    constexpr Point rot90() const {return Point(-y, x);}
    constexpr Point rot(long double ang) const {
        return Point(cos(ang) * x - sin(ang) * y, sin(ang) * x + cos(ang) * y);
    }
    
    // arithmetic operators
    constexpr Point operator - () const {return Point(-x, -y);}
    constexpr Point operator + (const Point &r) const {return Point(*this) += r;}
    constexpr Point operator - (const Point &r) const {return Point(*this) -= r;}
    constexpr Point operator * (const Point &r) const {return Point(*this) *= r;}
    constexpr Point operator / (const Point &r) const {return Point(*this) /= r;}
    constexpr Point operator * (DD r) const {return Point(*this) *= r;}
    constexpr Point operator / (DD r) const {return Point(*this) /= r;}
    constexpr Point& operator += (const Point &r) {
        x += r.x, y += r.y;
        return *this;
    }
    constexpr Point& operator -= (const Point &r) {
        x -= r.x, y -= r.y;
        return *this;
    }
    constexpr Point& operator *= (const Point &r) {
        DD tx = x, ty = y;
        x = tx * r.x - ty * r.y;
        y = tx * r.y + ty * r.x;
        return *this;
    }
    constexpr Point& operator /= (const Point &r) {
        return *this *= r.conj() / r.norm();
    }
    constexpr Point& operator *= (DD r) {
        x *= r, y *= r;
        return *this;
    }
    constexpr Point& operator /= (DD r) {
        x /= r, y /= r;
        return *this;
    }

    // friend functions
    friend ostream& operator << (ostream &s, const Point &p) {
        return s << '(' << p.x << ", " << p.y << ')';
    }
    friend constexpr Point conj(const Point &p) {return p.conj();}
    friend constexpr DD dot(const Point &p, const Point &q) {return p.dot(q);}
    friend constexpr DD cross(const Point &p, const Point &q) {return p.cross(q);}
    friend constexpr DD norm(const Point &p) {return p.norm();}
    friend constexpr long double abs(const Point &p) {return p.abs();}
    friend constexpr long double amp(const Point &p) {return p.amp();}
    friend constexpr bool eq(const Point &p, const Point &q) {return p.eq(q);}
    friend constexpr Point rot90(const Point &p) {return p.rot90();}
    friend constexpr Point rot(const Point &p, long long ang) {return p.rot(ang);}
};

// Line
struct Line : vector<Point> {
    Line(Point a = Point(0.0, 0.0), Point b = Point(0.0, 0.0)) {
        this->push_back(a);
        this->push_back(b);
    }
    friend ostream& operator << (ostream &s, const Line &l) {
        return s << '{' << l[0] << ", " << l[1] << '}';
    }
};


/*/////////////////////////////*/
// 直線の交差判定, 距離
/*/////////////////////////////*/

int ccw_for_dis(const Point &a, const Point &b, const Point &c) {
    if (cross(b-a, c-a) > EPS) return 1;
    if (cross(b-a, c-a) < -EPS) return -1;
    if (dot(b-a, c-a) < -EPS) return 2;
    if (norm(b-a) < norm(c-a) - EPS) return -2;
    return 0;
}
Point proj(const Point &p, const Line &l) {
    DD t = dot(p - l[0], l[1] - l[0]) / norm(l[1] - l[0]);
    return l[0] + (l[1] - l[0]) * t;
}
Point refl(const Point &p, const Line &l) {
    return p + (proj(p, l) - p) * 2;
}
bool is_inter_PL(const Point &p, const Line &l) {
    return (abs(p - proj(p, l)) < EPS);
}
bool is_inter_PS(const Point &p, const Line &s) {
    return (ccw_for_dis(s[0], s[1], p) == 0);
}
bool is_inter_LL(const Line &l, const Line &m) {
    return (abs(cross(l[1] - l[0], m[1] - m[0])) > EPS ||
            abs(cross(l[1] - l[0], m[0] - l[0])) < EPS);
}
bool is_inter_SS(const Line &s, const Line &t) {
    if (eq(s[0], s[1])) return is_inter_PS(s[0], t);
    if (eq(t[0], t[1])) return is_inter_PS(t[0], s);
    return (ccw_for_dis(s[0], s[1], t[0]) * ccw_for_dis(s[0], s[1], t[1]) <= 0 &&
            ccw_for_dis(t[0], t[1], s[0]) * ccw_for_dis(t[0], t[1], s[1]) <= 0);
}
DD distance_PL(const Point &p, const Line &l) {
    return abs(p - proj(p, l));
}
DD distance_PS(const Point &p, const Line &s) {
    Point h = proj(p, s);
    if (is_inter_PS(h, s)) return abs(p - h);
    return min(abs(p - s[0]), abs(p - s[1]));
}
DD distance_LL(const Line &l, const Line &m) {
    if (is_inter_LL(l, m)) return 0;
    else return distance_PL(m[0], l);
}
DD distance_SS(const Line &s, const Line &t) {
    if (is_inter_SS(s, t)) return 0;
    else return min(min(distance_PS(s[0], t), distance_PS(s[1], t)),
                    min(distance_PS(t[0], s), distance_PS(t[1], s)));
}

int main() {
    int N;
    cin >> N;
    vector<Line> segs(N);
    for (int i = 0; i < N; ++i) {
        int a, b, c, d;
        cin >> a >> b >> c >> d;
        segs[i] = Line(Point(a, b), Point(c, d));
    }
    
    // (x, y) を中心としたときの必要最小半径
    auto f = [&](DD x, DD y) -> DD {
        DD res = 0;
        for (const auto &seg : segs) {
            DD dis = distance_PS(Point(x, y), seg);
            res = max(res, dis);
        }
        return res;
    };
    
    // x を固定したときの f(x, y) の最大値を求める
    auto max_y_f = [&](DD x) -> DD {
        DD left = -1000, right = 3000;
        for (int iter = 0; iter < 100; ++iter) {
            DD m1 = (left * 2 + right) / 3;
            DD m2 = (left + right * 2) / 3;
            if (f(x, m1) > f(x, m2)) left = m1;
            else right = m2;
        }
        return f(x, right);
    };
    
    // f(x, y) の最大値を求める
    auto max_xy_f = [&]() -> DD {
        DD left = -1000, right = 3000;
        for (int iter = 0; iter < 100; ++iter) {
            DD m1 = (left * 2 + right) / 3;
            DD m2 = (left + right * 2) / 3;
            if (max_y_f(m1) > max_y_f(m2)) left = m1;
            else right = m2;
        }
        return max_y_f(right);
    };
    
    cout << fixed << setprecision(10) << max_xy_f() << endl;
}