BZOJ 2178 圆的面积并 ——Simpson积分

【题目分析】

    史上最良心样例，史上最难调样例。

    Simpson积分硬上。

    听说用long double 精度1e-10才能过。

    但是double+1e-6居然过了。

【代码】

#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
 
#include <map>
#include <set>
#include <queue>
#include <string>
#include <iostream>
#include <algorithm>
 
using namespace std;
 
#define maxn 1005
#define eps 1e-6
#define db double
#define ll long long
#define ldb long double
#define inf 0x3f3f3f3f
#define F(i,j,k) for (int i=j;i<=k;++i)
#define D(i,j,k) for (int i=j;i>=k;--i)
 
void Finout()
{
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
//    freopen("out.txt","w",stdout);
    #endif
}
 
int Getint()
{
    int x=0,f=1; char ch=getchar();
    while (ch<'0'||ch>'9') {if (ch=='-') f=-1; ch=getchar();}
    while (ch>='0'&&ch<='9') {x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}

int n,cnt=0,tag[maxn],tot;
struct Circle{int x,y,r;}c[maxn];
struct Segment{db l,r;}a[maxn];

bool cmp(Segment x,Segment y)
{return x.l==y.l?x.r<y.r:x.l<y.l;}

bool cmp1(Circle a,Circle b)
{return a.r<b.r;}

db get(db x)
{
	int cnt=0;
	F(i,1,n) if (fabs(c[i].x-x)<=c[i].r-eps)
	{
		db tmp=sqrt(c[i].r*c[i].r-(x-c[i].x)*(x-c[i].x));
		a[++cnt]=(Segment){c[i].y-tmp,c[i].y+tmp};
	}
	sort(a+1,a+cnt+1,cmp);
	db h=-inf,ans=0;
	F(i,1,cnt)
	{
		if (a[i].l>h) ans+=a[i].r-a[i].l,h=a[i].r;
		else if (a[i].r>h) ans+=a[i].r-h,h=a[i].r;
	}
	return ans;
}

db cal(db l,db r)
{
	db mid=(l+r)/2,fl=get(l),fr=get(r),fm=get(mid);
	return (r-l)/6*(fl+fr+4*fm);
}

db simpson(db l,db r)
{
	db mid=(l+r)/2,s1=cal(l,r),s2=cal(l,mid)+cal(mid,r);
	if (fabs(s2-s1)<=eps) return s2;
	else return simpson(l,mid)+simpson(mid,r); 
}

int main()
{
	Finout();n=Getint();
	F(i,1,n)
	{
		c[i].x=Getint();
		c[i].y=Getint();
		c[i].r=Getint();
	}
	sort(c+1,c+n+1,cmp1);
	F(i,1,n-1) F(j,i+1,n)
	{
		if ((c[i].x-c[j].x)*(c[i].x-c[j].x)+(c[i].y-c[j].y)*(c[i].y-c[j].y)<=(c[j].r-c[i].r)*(c[j].r-c[i].r))
		{
			tag[i]=1;
			break;
		}
	}
	F(i,1,n) if (!tag[i]) c[++tot]=c[i];
	n=tot;
	printf("%.3f\n",simpson(-2000.0,2000.0));
}