B. Forcefield

题意:平面上 n100n\le 100 个圆,半径可能不同,但保证互不重叠。你需要用一根皮筋把所有的圆全都包住,要求皮筋的周长最小,同时求出此时皮筋包住的区域的面积。

本场一道比较好写的计算几何。首先考虑两个圆 A,BA,B 的情况,一定是 A,BA,B 的一段弧再加上两条外公切线,就能把这两个圆包住,同时长度最小。

我们把这段皮筋分成两个部分:直线部分和弧线部分。直线部分就是外公切线,弧线就是圆弧。

所以我们把它扩展到多个圆的情况:我们先两两求出其外公切,然后对这些点我们求出凸包。这样,我们得到的凸包全都是外公切点,相邻两点之间的连线要么是圆的割线,要么是两个圆之间的外公切线。

现在假设我们的凸包是按顺时针,我们先算周长 CC。每次拿出相邻两点 a,ba,b,如果 a,ba,b 属于不同的圆,则他们构成外公切线,即 CC+abC\gets C+|ab|;否则找出 abab 这条割线对应的弧,如果其对应的圆心在凸包里面,则对应的是劣弧,否则对应的是优弧。

这里有一个 trick。令凸包上的点是顺时针排列的,那么考虑凸包上相邻两点 vi,vi+1v_i,v_{i+1} 和其对应的圆心 cc,考虑叉积 V=vic×vi+1cV=\vec{v_ic}\times\vec{v_{i+1}c}:如果 V0V\ge 0,那么圆心在凸包外;否则 V>0V\gt 0,则圆心在凸包内。

面积也是类似的。我们先计算凸包的面积,然后如果相邻两点 a,ba,b 属于同一个圆,如果圆心 cc 在凸包内部,则先减去 abcabc 的面积,再加上一段劣弧对应的扇形面积;如果在凸包外部,则先加上 abcabc 三角形的面积,再加上优弧对应的面积。

Code

注意几点:

  1. 精度问题:long doubleset_precision(20)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include <bits/stdc++.h>
#include <iomanip>
#define M_PI 3.141592653589793
#define double long double
using namespace std;
constexpr int N = 1000;
constexpr double eps = 1e-9;

struct point {
double x, y;
int id;
};
point operator-(point a, point b) { return {a.x - b.x, a.y - b.y, -1}; }
double scale(point a) { return sqrt(a.x * a.x + a.y * a.y); }
double dot(point a, point b) { return a.x * b.x + a.y * b.y; }
double cross(point a, point b) { return a.x * b.y - a.y * b.x; }

int sign(double l, double r) {
if (l - r > eps) return 1;
else if (r - l > eps) return -1;
else return 0;
}

struct circle {
point c;
double r;
int id;

point coord(double ang) {
double cs = cos(ang), ss = sin(ang);
return {c.x + r * cs, c.y + r * ss, id};
}
};

vector<point> p;
circle c[N];
int n;

tuple<vector<point>, vector<point>> tangent(circle A, circle B) {
vector<point> a, b;
if (A.r < B.r) swap(A, B), swap(a, b);

double dist = sqrt((A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y));
double rsum = A.r + B.r, rdif = A.r - B.r;
double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
double ang = acos(rdif / dist);
a.push_back(A.coord(base + ang)), b.push_back(B.coord(base + ang));
a.push_back(A.coord(base - ang)), b.push_back(B.coord(base - ang));

double ang1 = acos(rsum / dist);
a.push_back(A.coord(base + ang1)), b.push_back(B.coord(base + ang1 + M_PI));
a.push_back(A.coord(base - ang1)), b.push_back(B.coord(base - ang1 + M_PI));

return {a, b};
}

int main() {
#ifdef ONLINE_JUDGE
std::cin.tie(0)->sync_with_stdio(0);
#endif
std::cin >> n;
for (int i = 1; i <= n; i++) std::cin >> c[i].c.x >> c[i].c.y >> c[i].r, c[i].id = i;

for (int i = 1; i <= n; i++) {
for (int j = i + 1; j <= n; j++) {
auto [x, y] = tangent(c[i], c[j]);
for (auto v : x) p.push_back(v);
for (auto v : y) p.push_back(v);
}
}

// convex
auto v = [&] {
sort(p.begin(), p.end(), [&](point a, point b) {
if (sign(a.x, b.x) != 0) return a.x < b.x;
else return a.y < b.y;
});
auto crs = [&](point p, point a, point b) { return cross(a - p, b - p); };

vector<point> stack;
int top = -1, m = p.size();
for (int i = 0; i < m; i++) {
while (top >= 1 && sign(crs(stack[top - 1], stack[top], p[i]), 0) <= 0) {
stack.pop_back();
top--;
}
stack.push_back(p[i]), top++;
}

int t = top;
for (int i = m - 2; i >= 0; i--) {
while (top > t && sign(crs(stack[top - 1], stack[top], p[i]), 0) <= 0) {
stack.pop_back();
top--;
}
stack.push_back(p[i]);
top++;
}
return stack;
}();
std::reverse(v.begin(), v.end());

double C = 0;
double A = 0;
// compute A
for (int i = 0, m = v.size(); i + 1 < m; i++) A -= cross(v[i], v[i + 1]) * 0.5;

for (int i = 0, m = v.size(); i + 1 < m; i++) {
if (v[i].id != v[i + 1].id) {
C += scale(v[i] - v[i + 1]);
continue;
}
int cid = v[i].id;
double len = scale(v[i] - v[i + 1]);
double val = (c[cid].r * c[cid].r * 2 - len * len) / (c[cid].r * c[cid].r * 2);
double cosv = acos(min(1.0L, max(-1.0L, val)));

point v1 = v[i] - c[cid].c;
point v2 = v[i + 1] - c[cid].c;

A += cross(v[i] - c[cid].c, v[i + 1] - c[cid].c) * 0.5;

if (sign(cross(v[i] - c[cid].c, v[i + 1] - c[cid].c), 0) >= 0) {
cosv = 2 * M_PI - cosv;
C += c[cid].r * cosv;
A += cosv * c[cid].r * c[cid].r * 0.5;
} else C += c[cid].r * cosv, A += cosv * c[cid].r * c[cid].r * 0.5;
}
std::cout << setprecision(20) << fixed << C << " " << A << "\n";
}