#include <bits/stdc++.h>
using namespace std;
// #define double long double
const double eps = numeric_limits<double>::epsilon();
double _fabs(double x) {
if (x > -eps && x < eps) {
return 0;
} else {
return (x < 0) ? (-x) : x;
}
}
bool epse(double a, double b) {
return (_fabs(a - b) < eps) ? true : false;
}
struct Segment {
double x, y;
};
struct Coordinate {
double x, y;
double operator-(Coordinate const& a) const {
return sqrt(pow(this->x - a.x, 2) + pow(this->y - a.y, 2));
}
};
struct Circle {
Coordinate pos;
double r;
int relation(Circle x) {
double d = this->pos - x.pos;
if (epse(d, this->r + x.r)) {
return 0;
} else if (d > this->r + x.r) {
return 1;
} else {
return -1;
}
}
bool in(Coordinate p) {
double d = this->pos - p;
if (d < this->r) {
return true;
} else {
return false;
}
}
};
int n;
int cnt;
double ans;
vector<Circle> c;
vector<Circle> c_tmp;
bool cover[1005];
int st, ed;
vector<Segment> seg;
double f(double x) {
seg.clear();
for (int i = st; i <= ed; i++) {
if (x > c[i].pos.x + c[i].r || x < c[i].pos.x - c[i].r) {
continue;
}
double l = x - c[i].pos.x;
l = sqrt(c[i].r * c[i].r - l * l);
seg.push_back({c[i].pos.y - l, c[i].pos.y + l});
}
if (seg.empty()) {
return 0;
}
sort(seg.begin(), seg.end(), [](Segment a, Segment b) {
if (a.x == b.x) return a.y < b.y;
return a.x < b.y;
});
double ans = 0, last = seg[0].x;
for (int i = 0; i < seg.size(); ++i) {
if (seg[i].y > last) {
ans += seg[i].y - max(last, seg[i].x);
last = seg[i].y;
}
}
return ans;
}
double calc(double l, double r, double sl, double sr) {
double mid = (l + r) / 2.0;
return (r - l) * (sl + 4.0 * f(mid) + sr) / 6.0;
}
double simpson(double l, double r, double S, double sl, double sr) {
double mid = (l + r) / 2.0, s = f(mid);
double t1 = calc(l, mid, sl, s), t2 = calc(mid, r, s, sr);
if (_fabs(t1 + t2 - S) < eps) { return S; }
++cnt;
return simpson(l, mid, t1, sl, s) + simpson(mid, r, t2, s, sr);
}
int main() {
// freopen("SP8073.in", "r", stdin);
// freopen("SP8073.out", "w", stdout);
cin >> n;
for (int i = 1; i <= n; i++) {
double x, y, z;
cin >> x >> y >> z;
c_tmp.push_back({{x, y}, z});
}
sort(c_tmp.begin(), c_tmp.end(), [](Circle a, Circle b) {
return a.r > b.r;
});
for (auto i = 0; i < c_tmp.size(); i++) {
if (!cover[i]) {
c.push_back(c_tmp[i]);
}
for (int j = i + 1; j < c_tmp.size(); j++) {
if (c_tmp[i].relation(c_tmp[j]) == -1) {
cover[j] = true;
}
}
}
sort(c.begin(), c.end(), [](Circle a, Circle b) {
return a.pos.x - a.r < b.pos.x - b.r;
});
double l, r;
int i = 0, j;
while (i < c.size()) {
l = c[i].pos.x - c[i].r;
r = c[i].pos.x + c[i].r;
for (j = i + 1; j < c.size(); j++) {
if (c[j].pos.x - c[j].r > r) {
break;
}
r = max(r, c[j].pos.x + c[j].r);
}
st = i, ed = j - 1, i = j;
double mid = (l + r) / 2.0;
ans += simpson(l, r, f(mid), f(l), f(r));
}
printf("%lf\n", ans);
// fclose(stdin);
// fclose(stdout);
// system("pause");
return 0;
}
自适应辛普森法,SPOJ 上在 14 - 16 测试点之间 WA。