Bzoj2178:圆的面积并:Simpson积分

题目链接：2178:圆的面积并

Simpson积分裸上即可Ac，精度必须开到1e-13,1e-12都会wa

貌似有高大上的做法但是不会

Simpson真是慢死了QAQ

这篇博文用来学习Simpson很好

#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;
const double eps=1e-13;
const int maxn=1010;
const int inf=1e9+7;
int n;
struct point{double x,y;}p[maxn];
struct Circle{double x,y,r;}o[maxn],tmp[maxn];
bool mark[maxn];

bool cmp(const Circle &a,const Circle &b){
	return a.x-a.r<b.x-a.r||(a.x-a.r==b.x-b.r&&a.x+a.r<b.x+a.r);
}

bool cmp2(const point &a,const point &b){
	return a.x<b.x;
}

double qdis(double a,double b,double c,double d){
	return sqrt((a-c)*(a-c)+(b-d)*(b-d));
}

double getlen(double x,int id){
	if (o[id].r<fabs(o[id].x-x)) return -1;
	double tmp=o[id].r*o[id].r-(o[id].x-x)*(o[id].x-x);
	tmp=sqrt(tmp);
	return tmp;
}

double getseq(double x){
	double ret=0,last=-inf;
	int cnt=0;
	for (int i=1;i<=n;++i){
		double tmp=getlen(x,i);
		if (tmp>=0) p[++cnt].x=o[i].y-tmp,p[cnt].y=o[i].y+tmp;
	}
	sort(p+1,p+cnt+1,cmp2);
	for (int i=1;i<=cnt;++i){
        if (p[i].x>last) ret+=p[i].y-p[i].x,last=p[i].y;
        else if (p[i].y>last) ret+=p[i].y-last,last=p[i].y;
	}
	return ret;
}

double Simpson(double l,double r,double m,double al,double ar,double am){
	double t1=getseq((l+m)/2),t2=getseq((m+r)/2);
	double a1=(r-l)*(al+ar+4*am)/6,a2=(m-l)*(al+am+4*t1)/6,a3=(r-m)*(am+ar+4*t2)/6;
	if (fabs(a2+a3-a1)<eps) return a1;
	else return Simpson(l,m,(l+m)/2,al,am,t1)+Simpson(m,r,(m+r)/2,am,ar,t2);
}

int main(){
	scanf("%d",&n);
	double l=inf,r=-inf;
	for (int i=1;i<=n;++i){
		scanf("%lf%lf%lf",&o[i].x,&o[i].y,&o[i].r);
		l=min(l,o[i].x-o[i].r);
		r=max(r,o[i].x+o[i].r);
	}
	sort(o+1,o+n+1,cmp);
	for (int i=1;i<=n;++i){
		if (mark[i]) continue;
		for (int j=i+1;j<=n;++j){
			if (mark[j]) continue;
			if (qdis(o[i].x,o[i].y,o[j].x,o[j].y)+o[j].r<=o[i].r) mark[j]=1;
		}
	}
	int cnt=0;
	for (int i=1;i<=n;++i)
	    if (!mark[i]) tmp[++cnt]=o[i];
	for (int i=1;i<=cnt;++i) o[i]=tmp[i];
	n=cnt;
	printf("%.3lf",Simpson(l,r,(l+r)/2,0,0,getseq((l+r)/2)));
}