自适应Simpson积分

Posted 2020-10-24 探险家Mr.H

首先介绍一下Simpson积分的公式$$\int_{a}^{b}f(x)dx \approx \frac {b-a} {6} (f(a)+4f(\frac{a+b} {2})+f(b))$$

可以用这个求比较平滑的曲线的积分的近似...

但是直接用这个公式来做会误差很大

具体用法是：

用这个方法积$[l,mid]$，积$[mid,r]$，再积$[l,r]$判断是不是在误差范围内

如果不在就递归地细化$[l,mid]$和$[mid,r]$

比如说SDOI2017龙与地下城这道题

就用这种方法积了一个正态分布曲线

#include<bits/stdc++.h>

using namespace std;

inline int read()
{
    int x=0;char ch=getchar();
    while(ch<‘0‘ || ‘9‘<ch)ch=getchar();
    while(‘0‘<=ch && ch<=‘9‘)x=x*10+(ch^48),ch=getchar();
    return x;
}

typedef double db;

db x,y;

namespace Fast_Fouriet_Transform
{
    typedef complex<db> cp;
    const int M=800009;
    const db pi=acos(-1);

    cp a[M],b[M];
    db per,ans[M];
    int m,l,rev[M];

    inline void FFT(cp *a,int n,int f)
    {
        for(int i=0;i<n;i++)
            if(i<rev[i])
                swap(a[i],a[rev[i]]);
        for(int h=2;h<=n;h<<=1)
        {
            cp w(cos(2*pi/h),f*sin(2*pi/h));
            for(int j=0;j<n;j+=h)
            {
                cp wn(1,0),x,y;
                for(int k=j;k<j+(h>>1);k++)
                {
                    x=a[k],y=wn*a[k+(h>>1)];
                    a[k]=x+y;
                    a[k+(h>>1)]=x-y;
                    wn*=w;
                }
            }
        }
        if(f==-1)
            for(int i=0;i<n;i++)
                a[i].real()/=(db)n;
    }

    int mian()
    {
        per=1.0/(db)x;

        for(m=1,l=0;m<(y*x);m<<=1)l++;
        for(int i=0;i<m;i++)
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));

        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));

        a[0].real()=1.0;
        for(int i=0;i<x;i++)
            b[i].real()=per;

        FFT(a,m,1);FFT(b,m,1);
        for(int i=0;i<m;i++)
        {
            int cnt=y;
            while(cnt)
            {
                if(cnt&1)
                    a[i]=a[i]*b[i];
                b[i]=b[i]*b[i];
                cnt>>=1;
            }
        }
        FFT(a,m,-1);

        ans[0]=a[0].real();
        for(int i=1;i<m;i++)
            ans[i]=a[i].real()+ans[i-1];

        for(int i=1;i<=10;i++)
        {
            int a=read(),b=read();
            if(a==0)
                printf("%.8lf\n",ans[b]);
            else
                printf("%.8lf\n",ans[b]-ans[a-1]);
        }
        return 0;
    }
}

namespace simpsons
{
    const db pi2=sqrt(acos(-1)*2.0);
    const db eps=1e-12;
    db sigma,mu;
    inline db sqr(db x){return x*x;}
    inline db f(db x){return exp(-sqr(x)/2)/pi2;}
    inline bool dcmp(db x){return fabs(x)<15*eps;}
    inline void init()
    {
        sigma=sqrt((sqr(x)-1.0)/12.0);
        mu=(x-1.0)/2.0;
    }

    inline db simpson(db a,db b,db fl,db fmid,db fr)
    {
        return (b-a)*(fl+4.0*fmid+fr)/6.0;
    }

    inline db solve(db l,db r)
    {
        db mid=(l+r)/2.0,lm=(l+mid)/2.0,mr=(mid+r)/2.0;
        db fl=f(l),fmid=f(mid),fr=f(r),flm=f(lm),fmr=f(mr);

        db sipl=simpson(l,mid,fl,flm,fmid);
        db sipr=simpson(mid,r,fmid,fmr,fr);
        db sipm=simpson(l,r,fl,fmid,fr);

        if(dcmp(sipl+sipr-sipm))
            return sipl+sipr+(sipl+sipr-sipm)/15;
        else return solve(l,mid)+solve(mid,r);
    }

    int mian()
    {
        init();
        db base=solve(0,(x-1)*y);
        for(int i=1;i<=10;i++)
        {
            db a=(((db)read()-0.5)/y-mu)*sqrt(y)/sigma;
            db b=(((db)read()+0.5)/y-mu)*sqrt(y)/sigma;
            printf("%.8lf\n",solve(0,b)-solve(0,a));
        }
    }
}

int main()
{
    int T=read();
    while(T--)
    {
        scanf("%lf%lf",&x,&y);
        if(x*y<=400000)
            Fast_Fouriet_Transform::mian();
        else
            simpsons::mian();
    }

    return 0;
}

View Code

还有经典的NOI2005月下柠檬树

        L=min(L,C[i].c.x-C[i].r),R=max(R,C[i].c.x+C[i].r); 
//  printf("%lf\n%lf\n",L,R); 
    for (int i=1; i<=N; i++) 
        { 
            double d=C[i+1].c.x-C[i].c.x;  
            if (dcmp(d-fabs(C[i].r-C[i+1].r))<0) continue; //特判小圆被大圆覆盖的情况
            double sina=(C[i].r-C[i+1].r)/d,cosa=sqrt(1-sina*sina); 
            l[++num]=(Line){(Point){C[i].c.x+C[i].r*sina,C[i].r*cosa},(Point){C[i+1].c.x+C[i+1].r*sina,C[i+1].r*cosa}}; 
        } 
    printf("%.2lf\n",Simpson(L,R,Calc(L,R))); 
} 
int main() 
{ 
    scanf("%d%lf",&N,&alpha); 
    double h,r; 
    for (int i=1; i<=N+1; i++) 
        scanf("%lf",&h), 
        C[i]=(Circle){((Point){(h/tan(alpha))+C[i-1].c.x,0}),0}; 
    for (int i=1; i<=N; i++) 
        scanf("%lf",&r),C[i].r=r; 
//  for (int i=1; i<=N+1; i++) 
//      printf("%d     %.2lf     %.2lf\n",i,C[i].c.x,C[i].r); 
    Solve(); 
    return 0; 
}

View Code

感觉...算个黑科技吧

攒着以后可能会有点小用