bzoj 2178 圆的面积并【辛普森积分】

解题思路：

辛普森积分可以求过三点(l,f(l)),(mid,f(mid)),(r,f(r))，其中mid=(l+r)/2的抛物线在(l,r)一段的积分。其公式为：

∫rlf(x)=(r−l)6(f(l)+4f(mid)+f(r))

在这道题中，我们可以将f(i)取成直线x=i与所有圆的截线总长(重叠部分只算一次)，就可以算(l,r)一段所截的面积了。

由于不是标准的抛物线形式，所以我们可以递归求解，直到a
14a39
nsl+ansr和ans的误差在eps以内。这道题eps要取到1e-13。

还有个优化就是先去掉被其他圆包含的圆。

#include<cstdio>
#include<iostream>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
#include<ctime>
#include<vector>
#include<set>
#include<complex>
#define ll long long
using namespace std;

int getint()
{
int i=0,f=1;char c;
for(c=getchar();(c<'0'||c>'9')&&c!='-';c=getchar());
if(c=='-')f=-1,c=getchar();
for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
return i*f;
}

const int N=1005;
const double eps=1e-13,INF=1e9;
struct point
{
double x,y;
point(){}
point(double _x,double _y):
x(_x),y(_y){}
inline friend point operator - (const point &a,const point &b)
{return point(a.x-b.x,a.y-b.y); }
inline double dis() {return sqrt(x*x+y*y);}
inline friend bool operator <(const point &a,const point &b)
{
if(fabs(a.x-b.x)<eps)return a.y>b.y;
return a.x<b.x;
}
}l
;
struct circle
{
point O;double r;
inline friend bool operator < (const circle &a,const circle &b)
{return a.r<b.r;}
}C
;
int n;

void pre()
{
sort(C+1,C+n+1);
int tmp=0;
for(int i=1,j;i<=n;i++)
{
for(j=i+1;j<=n;j++)
if((C[i].O-C[j].O).dis()+C[i].r-C[j].r<eps)break;
if(j==n+1)C[++tmp]=C[i];
}
n=tmp;
}

double Cut(double x)
{
int cnt=0;double len;
for(int i=1;i<=n;i++)
if(abs(C[i].O.x-x)<C[i].r)
{
len=C[i].O.x-x;
len=sqrt(C[i].r*C[i].r-len*len);
l[++cnt]=point(C[i].O.y-len,C[i].O.y+len);
}
sort(l+1,l+cnt+1);
int top=0;
for(int i=1;i<=cnt;i++)
if(!top||l[i].y>l[top].y)l[++top]=l[i];
len=0;
for(int i=1;i<=top;i++)
if(i==1||l[i].x>l[i-1].y)len+=l[i].y-l[i].x;
else len+=l[i].y-l[i-1].y;
return len;
}

double Simpson(double l,double mid,double r,double fl,double fmid,double fr)
{
double lenl=Cut((l+mid)/2),lenr=Cut((mid+r)/2);
double ans=(r-l)*(fl+4*fmid+fr)/6,ansl=(mid-l)*(fl+4*lenl+fmid)/6,ansr=(r-mid)*(fmid+4*lenr+fr)/6;
if(abs(ans-ansl-ansr)<eps)return ans;
return Simpson(l,(l+mid)/2,mid,fl,lenl,fmid)+Simpson(mid,(mid+r)/2,r,fmid,lenr,fr);
}

int main()
{
//freopen("lx.in","r",stdin);
//freopen("lx.out","w",stdout);
n=getint();
double l=INF,r=-INF;
for(int i=1;i<=n;i++)
{
C[i].O.x=getint(),C[i].O.y=getint(),C[i].r=getint();
l=min(l,C[i].O.x-C[i].r),r=max(r,C[i].O.x+C[i].r);
}
pre();
printf("%0.3f",Simpson(l,(l+r)/2,r,0,Cut((l+r)/2),0));
return 0;
}