您的位置:首页 > 其它

[BZOJ2178]圆的面积并(辛普森积分)

2017-01-18 16:36 197 查看

题目描述

传送门

题解

这里的辛普森积分就不是算函数值了,而是算覆盖到的长度

有一个横坐标了之后O(n)枚举圆算弦长然后做一遍线段覆盖就行了

一遍过样例感觉好爽

然而被卡得狂T不止

主要有两个优化:

①调用simpson的时候多记录几个f值在下一次少算一些。能快好多

②把内含的圆去掉。能快好多好多

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define N 1005

const double eps=1e-13;
int dcmp(double x)
{
if (x<=eps&&x>=-eps) return 0;
return (x>0)?1:-1;
}
struct Circle
{
double x,y,r;
Circle(double X=0,double Y=0,double R=0)
{
x=X,y=Y,r=R;
}
bool operator < (const Circle &a) const
{
return r<a.r;
}
};
struct Segment
{
double l,r;
Segment(double L=0,double R=0)
{
l=L,r=R;
}
bool operator < (const Segment &a) const
{
return l<a.l;
}
};
int n;
double x,y,r,ans,Min=1e3,Max=-1e3;
Circle c
;
Segment a
;
bool flag
;

double F(double x)
{
int cnt=0;
for (int i=1;i<=n;++i)
{
double dis=fabs(c[i].x-x);
if (dcmp(dis-c[i].r)<0)
{
double len=sqrt(c[i].r*c[i].r-dis*dis);
a[++cnt]=Segment(c[i].y-len,c[i].y+len);
}
}
if (!cnt) return 0;
sort(a+1,a+cnt+1);
double l=a[1].l,r=a[1].r,ans=0;
for (int i=2;i<=cnt;++i)
if (dcmp(a[i].l-r)<=0) r=max(r,a[i].r);
else ans+=r-l,l=a[i].l,r=a[i].r;
ans+=r-l;
return ans;
}
double simpson(double l,double r,double now,double fl,double fr,double fmid)
{
double mid=(l+r)/2,flmid=F((l+mid)/2),frmid=F((mid+r)/2);
double p=(fl+fmid+flmid*4)*(mid-l)/6,q=(fmid+fr+frmid*4)*(r-mid)/6;
if (dcmp(now-p-q)==0) return now;
else return simpson(l,mid,p,fl,fmid,flmid)+simpson(mid,r,q,fmid,fr,frmid);
}

int main()
{
scanf("%d",&n);
for (int i=1;i<=n;++i)
{
scanf("%lf%lf%lf",&x,&y,&r);
c[i]=Circle(x,y,r);
Min=min(Min,x-r);Max=max(Max,x+r);
}
sort(c+1,c+n+1);
for (int i=1;i<=n;++i)
for (int j=i+1;j<=n;++j)
{
if (dcmp( sqrt( (c[i].x-c[j].x)*(c[i].x-c[j].x)+(c[i].y-c[j].y)*(c[i].y-c[j].y) ) + c[i].r - c[j].r )<=0)
{
flag[i]=true;
break;
}
}
int m=0;
for (int i=1;i<=n;++i)
if (!flag[i]) c[++m]=c[i];
n=m;
double fl=F(Min),fr=F(Max),fmid=F((Min+Max)/2);
ans=simpson(Min,Max,(fl+fr+fmid*4)*(Max-Min)/6,fl,fr,fmid);
printf("%.3lf\n",ans);
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: