BZOJ 2178 圆的面积并

Description

给出N个圆，求其面积并

Input

先给一个数字N ,N< = 1000 接下来是N行是圆的圆心，半径，其绝对值均为小于1000的整数

Output

面积并，保留三位小数

Sample Input

721

…

Sample Output

12707279.690

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

simpson积分~

因为直接积分会掉，所以我们以横坐标为基准把相互重叠的圆分成一组，每组分别积分，这样可以大大降低时间复杂度。

具体积分方法是：

定义cal(l,r)=(num(l)+num(mid)*4+num(r))*(r-l)/6，如果

cal(l,mid)+cal(mid,r)==cal(l,r)（即差值<eps），就返回cal(l,mid)+cal(mid+r)，否则继续向左右递归。

另外在输入之后要预处理去掉所有能被其它圆完全覆盖的圆。

注意如果把lx和rx这两个属性也加入到struct node里面排序的话会T掉，比分出来记录要慢10s左右，实在是神奇啊~

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
#define eps 1e-13

int n,tot,st,ed;
double ans,l,r,xl[1001],xr[1001];

struct node{
	double x,y,r;
}a[1001],c[1001];

struct node1{
	double l,r;
}p[1001];

int read()
{
	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<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}

bool han(node u,node v)
{
	return (u.x-v.x)*(u.x-v.x)+(u.y-v.y)*(u.y-v.y)>(u.r-v.r)*(u.r-v.r) ? 0:1;
}

bool operator < (node u,node v)
{
	return u.r>v.r;
}

bool cmp(node u,node v)
{
	return u.x-u.r<v.x-v.r;
}

bool operator < (node1 u,node1 v)
{
	return u.l==v.l ? u.r<v.r:u.l<v.l;
}

double num(double u)
{
	tot=0;double now,l,r;int j;
	for(int i=st;i<=ed;i++)
	  if(xl[i]<=u && xr[i]>=u)
	  {
	  	now=sqrt(c[i].r*c[i].r-(c[i].x-u)*(c[i].x-u));
		p[++tot].l=c[i].y-now;p[tot].r=c[i].y+now;
	  }
	sort(p+1,p+tot+1);now=0;
	for(int i=1;i<=tot;i++)
	{
		l=p[i].l;r=p[i].r;
		for(j=i+1;j<=tot;j++)
		  if(p[j].l>r) break;
		  else r=max(r,p[j].r);
		now+=r-l;i=j-1;
	}
	return now;
}

double cal(double l,double mid,double r,double fl,double fmid,double fr,double s)
{
	double l1=(l+mid)/2,r1=(mid+r)/2,fl1=num(l1),fr1=num(r1);
	double g1=(fl+fmid+4*fl1)*(mid-l)/6,g2=(fmid+fr+4*fr1)*(r-mid)/6;
	if(fabs(g1+g2-s)<eps) return g1+g2;
	return cal(l,l1,mid,fl,fl1,fmid,g1)+cal(mid,r1,r,fmid,fr1,fr,g2);
}

int main()
{
	n=read();
	for(int i=1;i<=n;i++)
	{
		a[i].x=read(),a[i].y=read(),a[i].r=read();
	}
	sort(a+1,a+n+1);
	for(int i=1;i<=n;i++)
	{
		c[++tot]=a[i];
		for(int j=1;j<tot;j++)
		  if(han(a[j],a[i]))
		  {
		  	tot--;break;
		  }
	}
	n=tot;int j;
	sort(c+1,c+n+1,cmp);
	for(int i=1;i<=n;i++) xl[i]=c[i].x-c[i].r,xr[i]=c[i].x+c[i].r;
	for(int i=1;i<=n;i++)
	{
		l=xl[i];r=xr[i];
		for(j=i+1;j<=n;j++)
		  if(xl[j]>r) break;
		  else r=max(r,xr[j]);
		st=i;ed=j-1;i=j-1;
		double mid=(l+r)/2,numl=num(l),nummid=num(mid),numr=num(r);
		ans+=cal(l,mid,r,numl,nummid,numr,(numl+numr+nummid*4)*(r-l)/6);
	}
	printf("%.3lf\n",ans);
	return 0;
}