关于Simpson积分

这两天写了圆交和圆并

圆交和圆并都有非常优美的O(n^2logn)算法，AekdyCoin有讲

但是像这种求面积的题还可以用Simpson积分法

简单的说就是将一段函数积分用二次函数积分拟合

一听这种搞法就知道是乱搞……但是很多时候比较有用……

嗯……像是求圆并的话，可以这样做

将x轴某一点上各个圆并所对应的长度视为函数值，用Simpson积分拟合

圆并在x轴上某一点对应的长度比较好算

如图，蓝色的那一段就是这一点对应的函数值，那么实际上就是每个圆和直线求一次交，得到的一段线段再做一次线段覆盖

看的出来复杂度是比较高的……每求一次函数值就是O(n)

具体的拟合怎么做呢

对于一段区间(l,r)我们首先计算f(l),f(r),f(mid)

然后(l,r)对应的拟合积分就是（f(l)+f(r)+4*f(mid)）*(r-l)/6

然后再分别计算(l,mid),(mid,r)的积分

如果误差<eps，认为这一段拟合正确，否则递归计算两个区间

又好想又好写~~

问题是，这样做真是乱搞的……

在做得时候要注意几个问题

一是最好确保积的这一段函数是连续的，否则极有可能出现4个峰分别在mid两边，你积出来得0的情况……

二是注意精度……如果输出要求精确到几位，你的eps最好比他高几个数量级……否则容易出错……而精度设高的代价为非常慢……

综上所述，Simpson积分法虽然不失为一个非常优秀的算法

但是由于他的种种问题，在可以想出正解或者有比他更好的算法的时候尽量少用……

优点就是适用范围广，好想好实现

如果没想出什么更好的算法就尽管用吧，总比没分好……

但是如果是ACM或者TC，CF之类的需要全部AC的慎用……或者eps精度高一点……

话说这样的考试一般不会出积分才能过的题吧……

SPOJ CIRU的代码 //此题为裸圆并

//Lib#include<cstdio>#include<cstring>#include<cstdlib>#include<cmath>#include<ctime>#include<iostream>#include<algorithm>#include<vector>#include<string>#include<queue>#include<stack>#include<set>#include<map>using namespace std;//Macro#definerep(i,a,b)for(int i=a,tt=b;i<=tt;++i)#definedrep(i,a,b)for(int i=a,tt=b;i>=tt;--i)#defineerep(i,e,x)for(int i=x;i;i=e[i].next)#defineirep(i,x)for(__typeof(x.begin()) i=x.begin();i!=x.end();i++)#defineread()(strtol(ipos,&ipos,10))#definesqr(x)((x)*(x))#definepbpush_back#define mpmake_pair#definePSsystem("pause");typedeflong longll;typedefpair<int,int>pii;const int oo=~0U>>1;const double inf=1e100;const double eps=1e-6;string name="", in=".in", out=".out";//Varstruct CIR{double x,y,r;void init(){scanf("%lf%lf%lf",&x,&y,&r);}}o[1008],tmp[1008];int n,st,ed;double xl[1008],xr[1008],ans;pair<double,double> seg[1008];bool cmp1(const CIR &a,const CIR &b){return a.r>b.r;}bool cmp2(const CIR &a,const CIR &b){return a.x-a.r<b.x-b.r||fabs(a.x-a.r-(b.x-b.r))<=eps&&a.x+a.r<b.x+b.r;}bool In(int i,int j){return sqr(tmp[i].x-o[j].x)+sqr(tmp[i].y-o[j].y)<=sqr(tmp[i].r-o[j].r);}void Init(){scanf("%d",&n);rep(i,1,n)o[i].init();//,s[i]=i;sort(o+1,o+1+n,cmp1);int k=0,j;rep(i,1,n){for(j=1;j<=k;j++)if(In(j,i))break;if(j>k)tmp[++k]=o[i];}n=k;rep(i,1,n)o[i]=tmp[i];}double f(double x){int tot=0,j;double ret=0,dis,l,r;rep(i,st,ed){if(x<=xl[i]||x>=xr[i])continue;dis=sqrt(o[i].r-sqr(x-o[i].x));seg[++tot]=mp(o[i].y-dis,o[i].y+dis);}sort(seg+1,seg+1+tot);rep(i,1,tot){l=seg[i].first;r=seg[i].second;for(j=i+1;j<=tot;j++)if(seg[j].first>r)break;else r=max(seg[j].second,r);ret+=r-l;i=j-1;}return ret;}double Calc(double s,double fl,double fr,double fmid){return (fl+fr+4*fmid)*s/6;}double Simpson(double l,double mid,double r,double fl,double fm,double fr,double tot){double m1=(l+mid)*0.5,m2=(mid+r)*0.5;double fm1=f(m1),fm2=f(m2);double g1=Calc(mid-l,fl,fm,fm1),g2=Calc(r-mid,fm,fr,fm2);if(fabs(tot-g1-g2)<=eps)return g1+g2;return Simpson(l,m1,mid,fl,fm1,fm,g1)+Simpson(mid,m2,r,fm,fm2,fr,g2);}void Work(){sort(o+1,o+1+n,cmp2);rep(i,1,n)xl[i]=o[i].x-o[i].r,xr[i]=o[i].x+o[i].r,o[i].r*=o[i].r;double l,r,mid,fl,fr,fm;int j;rep(i,1,n){l=xl[i],r=xr[i];for(j=i+1;j<=n;j++)if(xl[j]>r)break;else r=max(xr[j],r);mid=(l+r)*0.5;st=i;ed=j-1;i=j-1;fl=f(l);fm=f(mid);fr=f(r);ans+=Simpson(l,mid,r,fl,fm,fr,Calc(r-l,fl,fr,fm));}printf("%.3lf\n",ans);}int main(){Init();Work();return 0;}

P.S:

bzoj上有道和这个一模一样的圆并，但是eps得设到1e-13才能AC，可见Simpson积分十分不稳定啊

P.S2:

bzoj要关了……不知道root出于什么原因……谨以此篇写的非常烂的不知道想干什么的文章纪念我的bzoj账号……