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-;

  LDB ans=;

  struct data{

      LDB x,y,r,va;

  }a[],tmp[];

  int st,en,n,del[];

  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=;

      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+,tmp+cnt+,mycomp);

    LDB ret=;int po;

    for (int i=;i<=cnt;i=po+){

      po=i;LDB last=*tmp[i].y-tmp[i].va;

      while (po<cnt&&tmp[po+].va<last){

          po++;

          last=max(last,*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+*fmid+fr)/);

  }

  void simpson(LDB l,LDB r,LDB fl,LDB fmid,LDB fr){

      LDB mid=(l+r)/;

      LDB mid1=(l+mid)/,mid2=(mid+r)/;

      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=;i<=n;i=po+){

        po=i,last=a[i].x+a[i].r;

      while (po<n&&a[po+].x-a[po+].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)/),getf(last));

    }

  }

  int main(){

    scanf("%d",&n);

    for (int i=;i<=n;i++) scanf("%lf%lf%lf",&tmp[i].x,&tmp[i].y,&tmp[i].r);

    for (int i=;i<=n;i++)

      for (int j=;j<=n;j++)

        if (i!=j)

        if (dis(tmp[i],tmp[j])<tmp[j].r-tmp[i].r) {del[i]=;break;}

    int cnt=;

    for (int i=;i<=n;i++)

      if (!del[i]) a[++cnt]=tmp[i],a[cnt].va=a[cnt].x-a[cnt].r;

    n=cnt;

    sort(a+,a+n+,mycomp);

    work();

    printf("%.3lf\n",ans);

  }