自适应Simpson积分

首先介绍一下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
",ans[b]);
            else
                printf("%.8lf
",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
",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;
}

还有经典的NOI2005月下柠檬树

#include<iostream> 
#include<cstdio> 
#include<cstring> 
#include<algorithm> 
#include<cmath> 
using namespace std; 
#define MAXN 1010 
double alpha; 
int N,num; 
#define INF 1e12 
#define eps 1e-5 
struct Point 
{ 
    double x,y; 
    Point (double X=0,double Y=0) {x=X; y=Y;} 
}; 
struct Circle 
{ 
    double r; 
    Point c; 
    Circle(Point C=(Point){0,0},double R=0) {c=C; r=R;} 
}C[MAXN]; 
struct Line 
{ 
    Point s,t; 
    double k,b; 
    Line(Point S=(Point){0,0},Point T=(Point){0,0})  
        { 
            s=S,t=T; 
            if (s.x>t.x) swap(s,t); 
            k=(s.y-t.y)/(s.x-t.x); 
            b=s.y-k*s.x; 
        } 
    double f(double x) {return k*x+b;} 
}l[MAXN]; 
int dcmp(double x) {if (fabs(x)<eps) return 0; return x<0? -1:1;} 
double F(double x) 
{ 
    double re=0; 
    for (int i=1; i<=N; i++) //枚举圆是否与扫描线有交
        { 
            double d=fabs(x-C[i].c.x); 
            if (dcmp(d-C[i].r)>0) continue; 
            double len=2*sqrt(C[i].r*C[i].r-d*d); 
            re=max(re,len);  
        } 
    for (int i=1; i<=num; i++) //枚举公切线
        if (x>=l[i].s.x && x<=l[i].t.x) re=max(re,2*l[i].f(x)); 
    return re; 
} //利用扫描线去判断
double Calc(double l,double r) {double mid=(l+r)/2; return (F(l)+F(r)+F(mid)*4)*(r-l)/6;} 
double Simpson(double l,double r,double now) 
{ 
    double mid=(l+r)/2; 
    double x=Calc(l,mid),y=Calc(mid,r); 
    if (!dcmp(now-x-y)) return now; 
    else return Simpson(l,mid,x)+Simpson(mid,r,y); 
} 
void Solve() 
{ 
    double L=INF,R=-INF; 
    for (int i=1; i<=N+1; i++) 
        L=min(L,C[i].c.x-C[i].r),R=max(R,C[i].c.x+C[i].r); 
//  printf("%lf
%lf
",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
",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
",i,C[i].c.x,C[i].r); 
    Solve(); 
    return 0; 
}

感觉...算个黑科技吧

攒着以后可能会有点小用