AtCoder ABC 157F Yakiniku Optimization Problem

题目链接：https://atcoder.jp/contests/abc157/tasks/abc157_f

题目大意

　　平面上有$N$个点，第$i$个点的坐标为$(x_i,~y_i)$，其权重为$c_i$。

　　现在要在平面上选一点$(X,~Y)$，设$A_i~=~c_i imes sqrt{left(X - x_i ight)^2 + left(Y-y_i ight)^2}$。

　　定义$T$为$A$数组中第$K$小的数。

　　问如何选点，能使$T$最小，输出最小值。

分析

　　二分$T$。

　　可以把$N$个点看成是$N$个圆盘，第$i$个圆盘圆心为$(x_i,~y_i)$，半径为$frac{T}{c_i}$。

　　然后二分的标准为：是否存在一块区域，叠了大于等于$K$个圆盘。

　　实际代码中，可以用点来代替区域。

　　暴力枚举所有可能的点即可。

　　如果两圆盘相交，那么交点就是可能的点。

　　如果两圆盘内含，那么可以选择圆心为代表点。

　　时间复杂度为$O(log(frac{c_{max}~~~*~~(x_{max}~~~+~~y_{max}~~)}{EPS})~*~N^3)$

　　关于如何计算两圆交点，参考博客A和博客B

代码如下

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

/*-------------------Define Start-------------------*/
typedef bool BL;                        // 布尔类型
typedef char SB;                        // 有符号1字节,8位
typedef unsigned char UB;                // 无符号1字节,8位
typedef short SW;                        // 有符号短整型,16位
typedef unsigned short UW;                // 无符号短整型,16位
typedef long SDW;                        // 有符号整型,32位
typedef unsigned long UDW;               // 无符号整型,32位
typedef long long SLL;                    // 有符号长整型,64位
typedef unsigned long long ULL;            // 无符号长整型,64位
typedef char CH;                        // 单个字符
typedef float R32;                        // 单精度浮点数
typedef double R64;                        // 双精度浮点数

#define Rep(i, n) for (register SDW i = 0; i < (n); ++i)
#define For(i, s, t) for (register SDW i = (s); i <= (t); ++i)
#define rFor(i, t, s) for (register SDW i = (t); i >= (s); --i)
#define foreach(i, c) for (__typeof(c.begin()) i = c.begin(); i != c.end(); ++i)
#define ms0(a) memset(a,0,sizeof(a))
#define msI(a) memset(a,0x7f,sizeof(a))
#define LOWBIT(x) ((x)&(-x))

#define MP make_pair
#define PB push_back
#define ft first
#define sd second
#define ALL(x) x.begin(),x.end()

#define pr(x) cout << #x << " = " << x << "  "
#define prln(x) cout << #x << " = " << x << endl

const ULL mod = 1e9 + 7;                //常用模数(可根据题目需要修改)
const ULL inf = 0x7fffffff;                //用来表示无限大
const ULL infLL = 0x7fffffffffffffffLL;    //用来表示无限大
const R64 EPS = 1e-8;

// 重载<<操作符，用来打印vector 
template < typename T >
ostream& operator<< (ostream& out, vector< T > vec) {
    foreach(i, vec) {
        if(i != vec.begin()) {
            out << " ";
        }
        out << *i;
    }
    return out;
}

// 浮点数正负零判断
SDW sgn(R64 x) {
    if(fabs(x) < EPS) {
        return 0;
    }
    return x > 0 ? 1 : -1;
}

inline R64 pow2(R64 x) {
    return x * x;
}
/*-------------------Define End-------------------*/

template<typename T>
struct Point{
    SDW id;
    T X,Y;
    Point< T >(T x = 0,T y = 0) : X(x), Y(y) {}
    inline T set(T x, T y) { return X = x, Y = y; }
    inline Point<T> unitVec() { return Point< T >(X/abs(X), Y/abs(Y)); }
    inline T centerDist2() { return X*X + Y*Y; }
    inline double centerDist() { return sqrt(X*X + Y*Y); }
    
    inline bool operator== (const Point<T>& x) const { return X == x.X && Y == x.Y; }
    inline bool operator< (const Point<T>& x) const {
        if(Y == x.Y)return X < x.X;
        return Y < x.Y;
    }
    inline bool operator> (const Point<T>& x)const{return x < *this;}
    
    bool between(const Point<T> &x, const Point<T> &y) {
        if(x.X == y.X && X == x.X && Y > min(x.Y, y.Y) && Y < max(x.Y, y.Y)) return true;
        if(x.Y == y.Y && Y == x.Y && X > min(x.X, y.X) && X < max(x.X, y.X)) return true;
        return false;
    }
        
    inline Point<T> operator* (const T& k) { return Point<T>(k*X, k*Y); }
    inline Point<T> operator/ (const T& k) { return Point<T>(X/k, Y/k); }
    inline Point<T> operator+ (const Point<T>& x) { return Point<T>(X+x.X,Y+x.Y); }
    inline Point<T> operator- (const Point<T>& x) { return Point<T>(X-x.X,Y-x.Y); }
    
    T operator^(const Point<T> &x) const { return X*x.Y - Y*x.X; }
      T operator*(const Point<T> &x) const { return X*x.X + Y*x.Y; }
};

template<typename T>
istream &operator>> (istream &in, Point<T> &x) {
    in >> x.X >> x.Y;
    return in;
}

template<typename T>
ostream &operator<< (ostream &out, const Point<T> &x) {
    out << "(" << x.X << ", " << x.Y << ")";
    return out;
}

template<typename T>
double dist(Point<T> a,Point<T> b) {
  return sqrt((a-b)*(a-b));
} 

struct Circle{
    Point< R64 > p; // Circle center
    R64 c, r;
};
typedef Circle Disk;

istream& operator>> (istream& in, Disk &d) {
    cin >> d.p; 
    in >> d.c;
    return in;
}

/*
    设两圆圆心分别为A,B 
    两圆半径分别为Ra,Rb 
    两圆交点分别为C(x1,y1),D(x2,y2),x1 <= x2
    AB,CD交于点E(x0,y0)
    AB的斜率为Kab
    CD的斜率为Kcd
    过点C作X轴的垂线，过点E作Y轴垂线，两直线交于点F 
*/
BL getIntersectionOfTwoCircles(Circle &C1, Circle &C2, R64 &x1, R64 &y1, R64 &x2, R64 &y2) {
    Point< R64 > A = C1.p;
    Point< R64 > B = C2.p;
    R64 Ra = C1.r, Rb = C2.r;
    
    R64 AB = dist(A, B); // 圆心距 
    if(sgn(Ra + Rb - AB) < 0 || sgn(fabs(Ra - Rb) - AB) > 0) { // 两个圆没有交点的情况 
        return false;
    }
    
    R64 AE = (pow2(Ra) - pow2(Rb) + pow2(AB)) / (2 * AB);
    R64 CE = sqrt(pow2(Ra) - pow2(AE)); 
    R64 x0, y0;
    R64 Kab, Kcd;
    R64 EF, CF;
    
    if(sgn(A.X - B.X) == 0) { // 交点连线平行于X轴 
        x0 = A.X;
        y0 = A.Y + (B.Y - A.Y) * (AE / AB);
        
        x1 = x0 - CE;
        y1 = y0;
        x2 = x0 + CE;
        y2 = y0;
    } 
    else if(sgn(A.Y - B.Y) == 0) { // 交点连线平行于Y轴 
        x0 = A.X + (B.X - A.X) * (AE / AB);
        y0 = A.Y;
        
        x1 = x0;
        y1 = y0 - CE;
        x2 = x0;
        y2 = y0 + CE;
    } 
    else {
        Kab = (A.Y - B.Y) / (A.X - B.X);
        Kcd = -1 / Kab;
        x0 = A.X + (B.X - A.X) * (AE / AB);
        y0 = A.Y + Kab * (x0 - A.X);
        EF = sqrt(pow2(CE) / (1 + pow2(Kcd)));
        
        x1 = x0 - EF;
        y1 = y0 + Kcd * (x1 - x0);
        x2 = x0 + EF;
        y2 = y0 + Kcd * (x2 - x0);
    }
    
    return true;
}

const UDW maxN = 6e1 + 7;
SDW N, K;
Disk disks[maxN];
R64 ans;

void input(){
    cin >> N >> K;
    Rep(i, N) {
        cin >> disks[i];
    }
}

// 计数有多少disk覆盖点p 
BL checkPoint(Point< R64 > p) {
    R64 dis;
    SDW cnt = 0;
    
    Rep(i, N) {
        dis = dist(p, disks[i].p);
        
        if(sgn(disks[i].r - dis) >= 0) {
            ++cnt;
        }
    }
    return cnt >= K;
}

// 查找是否有一块区域是K个或K个以上Disk的交集 
BL findAns(R64 T) {
    Rep(i, N) {
        disks[i].r = T / disks[i].c;
    }
    
    // 内含情况，枚举圆心
    Rep(i, N) {
        if(checkPoint(disks[i].p)) {
            return true;
        }
    }
    
    // 相交情况，枚举交点 
    Rep(i, N) {
        For(j, i + 1, N - 1) {
            R64 x1, y1, x2, y2;
            
            if(!getIntersectionOfTwoCircles(disks[i], disks[j], x1, y1, x2, y2)) {
                continue;
            }
            
            if(checkPoint(Point< R64 >(x1, y1)) || checkPoint(Point< R64 >(x2, y2))) {
                return true;
            }
        }
    }
    
    return false; 
}

void solve(){
    R64 L = 0;
    R64 R = 2e5;
    
    while(sgn(L - R) < 0) {
        R64 mid = (L + R) / 2;
        
        if(findAns(mid)) {
            R = mid;
        }
        else {
            L = mid + 1e-6;
        }
    }
    
    ans = R;
}

void output(){
    printf("%f", ans);
}

int main() {
    input();
    solve();
    output();
    return 0;
}