simpson

使用二次函数拟合复杂的连续函数求积分

对于（l,r）拟合的积分为（r-l）*(f(l)+4*f((l+r)/2)+f(r))/6

___________________________

BZOJ2178

#include <cstdio>
#include <cmath>
#include <algorithm>
#define LDB double
using namespace std;

  const LDB eps=1e-13;
  LDB ans=0;

  struct data{
      LDB x,y,r,va;
  }a[1001],tmp[1001];
  
  int st,en,n,del[1001];

  int mycomp (const data&a,const data&b){
      return(a.va<b.va);
  }
  
  LDB dis(data a,data b){
      return(sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)));
  }
  
  LDB getf(LDB xx){
      int cnt=0;
      for (int i=st;i<=en;i++)
        if (fabs(a[i].x-xx)<a[i].r-eps){
          tmp[++cnt]=a[i];tmp[cnt].va=a[i].y-sqrt(a[i].r*a[i].r-(a[i].x-xx)*(a[i].x-xx));    
      }
    sort(tmp+1,tmp+cnt+1,mycomp);
    
    LDB ret=0;int po;
    for (int i=1;i<=cnt;i=po+1){
      po=i;LDB last=2*tmp[i].y-tmp[i].va;
      while (po<cnt&&tmp[po+1].va<last){
          po++;
          last=max(last,2*tmp[po].y-tmp[po].va);
      }
      ret+=(last-tmp[i].va);    
    }
    return(ret);  
  }
  
  LDB cal(LDB l,LDB r,LDB fl,LDB fmid,LDB fr){
      return((r-l)*(fl+4*fmid+fr)/6);
  }
  
  void simpson(LDB l,LDB r,LDB fl,LDB fmid,LDB fr){
      LDB mid=(l+r)/2;
      LDB mid1=(l+mid)/2,mid2=(mid+r)/2;
      LDB f1=getf(mid1),f2=getf(mid2);
      LDB num1=cal(l,r,fl,fmid,fr);
      LDB num2=cal(l,mid,fl,f1,fmid)+cal(mid,r,fmid,f2,fr);
      if (fabs(num1-num2)<eps) {ans+=num1;return;}
      simpson(l,mid,fl,f1,fmid);simpson(mid,r,fmid,f2,fr);
  }
  
  void work(){
      int po;LDB last;
      for (int i=1;i<=n;i=po+1){
        po=i,last=a[i].x+a[i].r;
      while (po<n&&a[po+1].x-a[po+1].r<last+eps){
          po++;
          last=max(last,a[po].x+a[po].r);
      }
      st=i;en=po;
      simpson(a[i].x-a[i].r,last,getf(a[i].x-a[i].r),getf((a[i].x-a[i].r+last)/2),getf(last));    
    }
  }

  int main(){      
    scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%lf%lf%lf",&tmp[i].x,&tmp[i].y,&tmp[i].r);
    for (int i=1;i<=n;i++)
      for (int j=1;j<=n;j++)
        if (i!=j)
        if (dis(tmp[i],tmp[j])<tmp[j].r-tmp[i].r) {del[i]=1;break;}
    int cnt=0;
    for (int i=1;i<=n;i++)
      if (!del[i]) a[++cnt]=tmp[i],a[cnt].va=a[cnt].x-a[cnt].r;
    n=cnt;
    sort(a+1,a+n+1,mycomp);
    
    work();
    printf("%.3lf\n",ans);
  }