当数据较离散时,会导致错误;(及解决方法)
原题Code:
/*
自适应辛普森积分法:
令f(x0):直线X = x0与所有圆的长度并
将f(X)在[-INF ~ INF]上面积分即为:所有圆的面积并
//注: 当区间比较离散时,运气不好的话(感觉第一次的时候很有可能)l和r的x轴都没与圆有交集,
mid也与圆没有交集,那直接就等于0<eps了,直接G了。
这里可以打个补丁: 可以先将所有圆投影到x轴上,然后做个区间合并,
然后对每个区间分别做辛普森。
*/
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
typedef double db;
const int N = 1010;
const db eps = 1e-8;
int n;
struct Circle
{
PDD o;
db r;
}circle[N];
PDD q[2010];
int dcmp(db x, db y)
{
if(fabs(x - y) < eps) return 0;
if(x < y) return -1;
return 1;
}
db f(db x)
{
int cnt = 0;
for(int i = 0; i < n; i ++)
{
db a = circle[i].o.x, b = circle[i].o.y, R = circle[i].r;
db X = fabs(a - x);
if(dcmp(X, R) < 0)
{
db Y = sqrt(R * R - X * X);
q[cnt ++] = {b - Y, b + Y};
}
}
if(!cnt) return 0;
db res = 0;
sort(q, q + cnt);
db st = q[0].x, ed = q[0].y;
for(int i = 1; i < cnt; i ++)
if(q[i].x <= ed) ed = max(ed, q[i].y);
else
{
res += ed - st;
st = q[i].x, ed = q[i].y;
}
return res + ed - st;
}
db simpson(db l, db r)
{
db mid = (l + r) / 2;
return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}
db asr(db l, db r, db s)
{
db mid = (l + r) / 2;
db left = simpson(l, mid), right = simpson(mid, r);
if(fabs(left + right - s) < eps) return right + left;
return asr(l, mid, left) + asr(mid, r, right);
}
int main()
{
scanf("%d", &n);
db l = 2000, r = -2000;
for(int i = 0; i < n; i ++)
{
scanf("%lf %lf %lf", &circle[i].o.x, &circle[i].o.y, &circle[i].r);
l = min(l, circle[i].o.x - circle[i].r), r = max(r, circle[i].o.x + circle[i].r);
}
printf("%.3lf\n", asr(l - 100, r + 100, simpson(l, r)));
return 0;
}
数据较为离散时:
这个数据就会出错:
2
53 -154 2.2203502244049247
351 -459 0
解决方法:
当区间比较离散时,运气不好的话(感觉第一次的时候很有可能)l和r的x轴都没与圆有交集,
mid也与圆没有交集,那直接就等于0<eps了,直接G了。
这里可以打个补丁: 可以先将所有圆投影到x轴上,然后做个区间合并,
然后对每个区间分别做辛普森。
改进Code:(但是$eps>1e-8$会wa最后一个点,有0.001误差;$eps <= 1e-9$会T,呜呜呜)
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<long double, long double> PDD;
typedef long double db;
const int N = 1010;
const db eps = 1e-8;
int n;
struct Circle
{
PDD o;
db r;
}circle[N];
PDD q[2010];
int dcmp(db x, db y)
{
if(fabs(x - y) < eps) return 0;
if(x < y) return -1;
return 1;
}
db f(db x)
{
int cnt = 0;
for(int i = 0; i < n; i ++)
{
db a = circle[i].o.x, b = circle[i].o.y, R = circle[i].r;
db X = fabs(a - x);
if(dcmp(X, R) < 0)
{
db Y = sqrt(R * R - X * X);
q[cnt ++] = {b - Y, b + Y};
}
}
if(!cnt) return 0;
db res = 0;
sort(q, q + cnt);
db st = q[0].x, ed = q[0].y;
for(int i = 1; i < cnt; i ++)
if(q[i].x <= ed) ed = max(ed, q[i].y);
else
{
res += ed - st;
st = q[i].x, ed = q[i].y;
}
return res + ed - st;
}
db simpson(db l, db r)
{
db mid = (l + r) / 2;
return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}
db asr(db l, db r, db s)
{
db mid = (l + r) / 2;
db left = simpson(l, mid), right = simpson(mid, r);
if(fabs(left + right - s) < eps) return right + left;
return asr(l, mid, left) + asr(mid, r, right);
}
int main()
{
scanf("%d", &n);
db l = 2000, r = -2000;
PDD p[2010];
int cnt = 0;
for(int i = 0; i < n; i ++)
{
scanf("%lf %lf %lf", &circle[i].o.x, &circle[i].o.y, &circle[i].r);
p[cnt ++] = {circle[i].o.x - circle[i].r, circle[i].o.x + circle[i].r};
}
db res = 0;
sort(p, p + cnt);
db st = p[0].x, ed = p[0].y;
for(int i = 1; i < cnt; i ++)
if(ed >= p[i].x) ed = max(ed, p[i].y);
else
{
res += asr(st, ed, simpson(st, ed));
st = q[i].x, ed = q[i].y;
}
if(st != ed) res += asr(st, ed, simpson(st, ed));
printf("%.3lf\n", res);
return 0;
}