[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-13using 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;elsereturn 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 secondusing 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);}