题目链接：https://atcoder.jp/contests/abc157/tasks/abc157_f
题目大意
平面上有$N$个点，第$i$个点的坐标为$(x_i,~y_i)$，其权重为$c_i$。
分析
关于如何计算两圆交点，参考博客A和博客B代码如下
1 #include <bits/stdc++.h>
2 using namespace std;
3
4 /*-------------------Define Start-------------------*/
5 typedef bool BL; // 布尔类型
6 typedef char SB; // 有符号1字节,8位
7 typedef unsigned char UB; // 无符号1字节,8位
8 typedef short SW; // 有符号短整型,16位
9 typedef unsigned short UW; // 无符号短整型,16位
10 typedef long SDW; // 有符号整型,32位
11 typedef unsigned long UDW; // 无符号整型,32位
12 typedef long long SLL; // 有符号长整型,64位
13 typedef unsigned long long ULL; // 无符号长整型,64位
14 typedef char CH; // 单个字符
15 typedef float R32; // 单精度浮点数
16 typedef double R64; // 双精度浮点数
17
18 #define Rep(i, n) for (register SDW i = 0; i < (n); ++i)
19 #define For(i, s, t) for (register SDW i = (s); i <= (t); ++i)
20 #define rFor(i, t, s) for (register SDW i = (t); i >= (s); --i)
21 #define foreach(i, c) for (__typeof(c.begin()) i = c.begin(); i != c.end(); ++i)
22 #define ms0(a) memset(a,0,sizeof(a))
23 #define msI(a) memset(a,0x7f,sizeof(a))
24 #define LOWBIT(x) ((x)&(-x))
25
26 #define MP make_pair
27 #define PB push_back
28 #define ft first
29 #define sd second
30 #define ALL(x) x.begin(),x.end()
31
32 #define pr(x) cout << #x << " = " << x << " "
33 #define prln(x) cout << #x << " = " << x << endl
34
35 const ULL mod = 1e9 + 7; //常用模数(可根据题目需要修改)
36 const ULL inf = 0x7fffffff; //用来表示无限大
37 const ULL infLL = 0x7fffffffffffffffLL; //用来表示无限大
38 const R64 EPS = 1e-8;
39
40 // 重载<<操作符，用来打印vector
41 template < typename T >
42 ostream& operator<< (ostream& out, vector< T > vec) {
43 foreach(i, vec) {
44 if(i != vec.begin()) {
45 out << " ";
46 }
47 out << *i;
48 }
49 return out;
50 }
51
52 // 浮点数正负零判断
53 SDW sgn(R64 x) {
54 if(fabs(x) < EPS) {
55 return 0;
56 }
57 return x > 0 ? 1 : -1;
58 }
59
60 inline R64 pow2(R64 x) {
61 return x * x;
62 }
63 /*-------------------Define End-------------------*/
64
65 template<typename T>
66 struct Point{
67 SDW id;
68 T X,Y;
69 Point< T >(T x = 0,T y = 0) : X(x), Y(y) {}
70 inline T set(T x, T y) { return X = x, Y = y; }
71 inline Point<T> unitVec() { return Point< T >(X/abs(X), Y/abs(Y)); }
72 inline T centerDist2() { return X*X + Y*Y; }
73 inline double centerDist() { return sqrt(X*X + Y*Y); }
74
75 inline bool operator== (const Point<T>& x) const { return X == x.X && Y == x.Y; }
76 inline bool operator< (const Point<T>& x) const {
77 if(Y == x.Y)return X < x.X;
78 return Y < x.Y;
79 }
80 inline bool operator> (const Point<T>& x)const{return x < *this;}
81
82 bool between(const Point<T> &x, const Point<T> &y) {
83 if(x.X == y.X && X == x.X && Y > min(x.Y, y.Y) && Y < max(x.Y, y.Y)) return true;
84 if(x.Y == y.Y && Y == x.Y && X > min(x.X, y.X) && X < max(x.X, y.X)) return true;
85 return false;
86 }
87
88 inline Point<T> operator* (const T& k) { return Point<T>(k*X, k*Y); }
89 inline Point<T> operator/ (const T& k) { return Point<T>(X/k, Y/k); }
90 inline Point<T> operator+ (const Point<T>& x) { return Point<T>(X+x.X,Y+x.Y); }
91 inline Point<T> operator- (const Point<T>& x) { return Point<T>(X-x.X,Y-x.Y); }
92
93 T operator^(const Point<T> &x) const { return X*x.Y - Y*x.X; }
94 T operator*(const Point<T> &x) const { return X*x.X + Y*x.Y; }
95 };
96
97 template<typename T>
98 istream &operator>> (istream &in, Point<T> &x) {
99 in >> x.X >> x.Y;
100 return in;
101 }
102
103 template<typename T>
104 ostream &operator<< (ostream &out, const Point<T> &x) {
105 out << "(" << x.X << ", " << x.Y << ")";
106 return out;
107 }
108
109 template<typename T>
110 double dist(Point<T> a,Point<T> b) {
111 return sqrt((a-b)*(a-b));
112 }
113
114 struct Circle{
115 Point< R64 > p; // Circle center
116 R64 c, r;
117 };
118 typedef Circle Disk;
119
120 istream& operator>> (istream& in, Disk &d) {
121 cin >> d.p;
122 in >> d.c;
123 return in;
124 }
125
126 /*
127 设两圆圆心分别为A,B
128 两圆半径分别为Ra,Rb
129 两圆交点分别为C(x1,y1),D(x2,y2),x1 <= x2
130 AB,CD交于点E(x0,y0)
131 AB的斜率为Kab
132 CD的斜率为Kcd
133 过点C作X轴的垂线，过点E作Y轴垂线，两直线交于点F
134 */
135 BL getIntersectionOfTwoCircles(Circle &C1, Circle &C2, R64 &x1, R64 &y1, R64 &x2, R64 &y2) {
136 Point< R64 > A = C1.p;
137 Point< R64 > B = C2.p;
138 R64 Ra = C1.r, Rb = C2.r;
139
140 R64 AB = dist(A, B); // 圆心距
141 if(sgn(Ra + Rb - AB) < 0 || sgn(fabs(Ra - Rb) - AB) > 0) { // 两个圆没有交点的情况
142 return false;
143 }
144
145 R64 AE = (pow2(Ra) - pow2(Rb) + pow2(AB)) / (2 * AB);
146 R64 CE = sqrt(pow2(Ra) - pow2(AE));
147 R64 x0, y0;
148 R64 Kab, Kcd;
149 R64 EF, CF;
150
151 if(sgn(A.X - B.X) == 0) { // 交点连线平行于X轴
152 x0 = A.X;
153 y0 = A.Y + (B.Y - A.Y) * (AE / AB);
154
155 x1 = x0 - CE;
156 y1 = y0;
157 x2 = x0 + CE;
158 y2 = y0;
159 }
160 else if(sgn(A.Y - B.Y) == 0) { // 交点连线平行于Y轴
161 x0 = A.X + (B.X - A.X) * (AE / AB);
162 y0 = A.Y;
163
164 x1 = x0;
165 y1 = y0 - CE;
166 x2 = x0;
167 y2 = y0 + CE;
168 }
169 else {
170 Kab = (A.Y - B.Y) / (A.X - B.X);
171 Kcd = -1 / Kab;
172 x0 = A.X + (B.X - A.X) * (AE / AB);
173 y0 = A.Y + Kab * (x0 - A.X);
174 EF = sqrt(pow2(CE) / (1 + pow2(Kcd)));
175
176 x1 = x0 - EF;
177 y1 = y0 + Kcd * (x1 - x0);
178 x2 = x0 + EF;
179 y2 = y0 + Kcd * (x2 - x0);
180 }
181
182 return true;
183 }
184
185 const UDW maxN = 6e1 + 7;
186 SDW N, K;
187 Disk disks[maxN];
188 R64 ans;
189
190 void input(){
191 cin >> N >> K;
192 Rep(i, N) {
193 cin >> disks[i];
194 }
195 }
196
197 // 计数有多少disk覆盖点p
198 BL checkPoint(Point< R64 > p) {
199 R64 dis;
200 SDW cnt = 0;
201
202 Rep(i, N) {
203 dis = dist(p, disks[i].p);
204
205 if(sgn(disks[i].r - dis) >= 0) {
206 ++cnt;
207 }
208 }
209 return cnt >= K;
210 }
211
212 // 查找是否有一块区域是K个或K个以上Disk的交集
213 BL findAns(R64 T) {
214 Rep(i, N) {
215 disks[i].r = T / disks[i].c;
216 }
217
218 // 内含情况，枚举圆心
219 Rep(i, N) {
220 if(checkPoint(disks[i].p)) {
221 return true;
222 }
223 }
224
225 // 相交情况，枚举交点
226 Rep(i, N) {
227 For(j, i + 1, N - 1) {
228 R64 x1, y1, x2, y2;
229
230 if(!getIntersectionOfTwoCircles(disks[i], disks[j], x1, y1, x2, y2)) {
231 continue;
232 }
233
234 if(checkPoint(Point< R64 >(x1, y1)) || checkPoint(Point< R64 >(x2, y2))) {
235 return true;
236 }
237 }
238 }
239
240 return false;
241 }
242
243 void solve(){
244 R64 L = 0;
245 R64 R = 2e5;
246
247 while(sgn(L - R) < 0) {
248 R64 mid = (L + R) / 2;
249
250 if(findAns(mid)) {
251 R = mid;
252 }
253 else {
254 L = mid + 1e-6;
255 }
256 }
257
258 ans = R;
259 }
260
261 void output(){
262 printf("%f", ans);
263 }
264
265 int main() {
266 input();
267 solve();
268 output();
269 return 0;
270 }
View Code