[Simpson积分 || 圆的离散化 几何] BZOJ 2178 圆的面积并

不知道正解为何物：http://mojijs.com/2015/09/206783/index.html

直接套用Simpson 函数值是一些线段的并

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#define PI acos(-1.0)
#define eps 1e-13
using namespace std;
typedef pair<double,double> abcd;

inline char nc()
{
	static char buf[100000],*p1=buf,*p2=buf;
	if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
	return *p1++;
}

inline void read(int &x){
	char c=nc(),b=1;
	for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
	for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

struct Point{
	int x,y;
	void read(){
		::read(x); ::read(y);
	}
	friend double dist(Point a,Point b){
        return sqrt((double)(a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
    }
};

struct Circle{
	Point p;
	int r;
	void read(){
		p.read(); ::read(r);
	}
    abcd f(double x)
    {
        if(r<=fabs(p.x-x)) return abcd(0,0);
        double t=r*r-(p.x-x)*(p.x-x);
        t=sqrt(t);
        return abcd(p.y-t,p.y+t);
    }
}C[1005];

int n;
int del[1005];
abcd L[1005];

inline double F(double x){
	double ret=0,last=-1e130;
	int cnt=0;
	for (int i=1;i<=n;i++)
	{
		L[++cnt]=C[i].f(x);
		if (L[cnt]==abcd(0,0))
			cnt--;
	}
	sort(L+1,L+cnt+1);
	for (int i=1;i<=cnt;i++)
	{
    	if(L[i].first>last)
            ret+=L[i].second-L[i].first,last=L[i].second;
        else if(L[i].second>last)
            ret+=L[i].second-last,last=L[i].second;
    }
    return ret;
}

inline double Simpson(double l,double r,double mid,double fl,double fr,double fm){
	double flm=F((l+mid)/2),fmr=F((mid+r)/2);
	double ret=(fl+4*fm+fr)*(r-l)/6,lret=(fl+4*flm+fm)*(mid-l)/6,rret=(fm+4*fmr+fr)*(r-mid)/6;
	if (fabs(ret-lret-rret)<eps)
		return ret;
	else
		return Simpson(l,mid,(l+mid)/2,fl,fm,flm)+Simpson(mid,r,(mid+r)/2,fm,fr,fmr);
}

int main()
{
	double l=1e130,r=-1e130;
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	read(n);
	for (int i=1;i<=n;i++)
	{
		C[i].read();
		l=min(l,(double)C[i].p.x-C[i].r);
		r=max(r,(double)C[i].p.x+C[i].r);
	}
	for (int i=1;i<=n;i++)
	{
		int flag=0;
		for (int j=1;j<=n && !flag;j++)
			if (i!=j && C[j].r>C[i].r && dist(C[i].p,C[j].p)<=C[j].r-C[i].r)
				flag=1;
		if (flag)
			del[i]=1;
	}
	int pnt=0;
	for (int i=1;i<=n;i++)
		if (!del[i])
			C[++pnt]=C[i];
	n=pnt;
	printf("%.3lf\n",Simpson(l,r,(l+r)/2,0,0,F((l+r)/2)));
	return 0;
}

UPD：

正解！

http://wenku.baidu.com/link?url=WfXnvfrpHsg3g4ViWDVWvCNuEvBJsemEGmor808WeRY3yHifiojI6o2zMfkU9SfYf0niYcPNa3oHiH2S6KvEjVtCqV607lu-deYuWSL2PqK

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#define L first
#define R second
using namespace std;
typedef long double ld;
typedef pair<ld,ld> abcd;

const int N=1005;

inline ld sqr(ld x){ return x*x; }
const ld PI=acos(-1.0);
const ld eps=1e-15;
inline int dcmp(ld a,ld b){
	if (fabs(a-b)<eps) return 0;
	return a<b?-1:1;
}

inline ld K(ld x,ld y){
	x+=eps;
	if (x<0) return atan(y/x)+PI;
	if (y<0) return atan(y/x)+PI*2;
	return atan(y/x); 
}

struct Circle{
	ld x,y,r;
	friend ld dist(Circle A,Circle B){
		return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));
	}
	void read(){
		double _x,_y,_r; scanf("%lf%lf%lf",&_x,&_y,&_r); x=_x,y=_y,r=_r;
	}
	bool Cover(Circle B){
		return dcmp(B.r+dist(*this,B),r)<=0;
	}
}C[N];

int n; ld ans;
int cover[N];
abcd S[N*N]; int cnt;
int l,r;

int main()
{
	ld ix,iy,dd,tmp; ld x_1,x_2,y_1,y_2;
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	scanf("%d",&n);
	for (int i=1;i<=n;i++) C[i].read();
	for (int i=1;i<=n;i++)
	{
		cnt=0;
		for (int j=1;j<=n && !cover[i];j++){
			if (i==j || cover[j]) continue;
			if (C[j].Cover(C[i])) { cover[i]=1; continue; }
			if (C[i].Cover(C[j])) continue;
			if (dcmp(dd=dist(C[i],C[j]),C[i].r+C[j].r)<0){
				++cnt;
				tmp=K(C[j].x-C[i].x,C[j].y-C[i].y);  
				ix=(sqr(C[i].r)-sqr(C[j].r)+sqr(dd))/2/dd;  
				iy=sqrt(sqr(C[i].r)-ix*ix);
				S[cnt].L=tmp-K(ix,iy);  
				S[cnt].R=tmp+K(ix,iy);  
				if (S[cnt].L<0)
			  	S[cnt].L+=2*PI,S[cnt].R+=2*PI; 
			}
		}
		if (cnt==0 && !cover[i]) ans+=sqr(C[i].r)*PI;
		if (cover[i]) continue;//
		sort(S+1,S+cnt+1); 
		r=0,l=1;
		for (int j=1;j<=cnt;j++) 
			if (r==0 || dcmp(S[j].L,S[r].R)>0)
	  			S[++r]=S[j];
			else if (dcmp(S[r].R,S[j].R)<0)
				S[r].R=S[j].R;  
		for (;l<=r && dcmp(S[r].R-2*PI,S[l].L)>0;l++)
			if (dcmp(S[r].R-2*PI,S[l].R)<0) 
				S[r].R=S[l].R+2*PI;
		for (int j=l,k;j<=r;j++){
			j==r?k=l:k=j+1;
			ix=S[k].L-S[j].R; if (ix<0) ix+=2*PI;  
			ans+=sqr(C[i].r)*(ix-sin(ix))/2;  
   			x_1=C[i].x+C[i].r*cos(S[j].R);  
			y_1=C[i].y+C[i].r*sin(S[j].R);  
			x_2=C[i].x+C[i].r*cos(S[k].L);  
			y_2=C[i].y+C[i].r*sin(S[k].L);  
			ans+=(x_1*y_2-x_2*y_1)/2;  
		}
	}
	printf("%.3lf\n",(double)ans);
}