不知道正解为何物: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:
正解!
#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);
}