2014-team3/code/adaptive_simpson

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

{{{
typedef double flt;
flt adaptive(flt l,flt r,flt m,flt eps,flt s,flt fl,flt fr,flt fm,int lim){
    flt h=(r-l)/12;
    flt u=ldexp(l+m,-1),v=ldexp(m+r,-1);
    flt fu=f(u),fv=f(v);
    flt sl=(fl+fm+ldexp(fu,2))*h,
        sr=(fm+fr+ldexp(fv,2))*h,
        s2=sl+sr;
    if(lim<=0 || fabs(s2-s)<= eps)
        return s2+(s2-s)/15;
    eps=ldexp(eps,-1);--lim;
    return  adaptive(l,m,u,eps,sl,fl,fm,fu,lim)
        +   adaptive(m,r,v,eps,sr,fm,fr,fv,lim);
}
flt simpson(flt l,flt r,flt eps=1e-6,int maxrec=10){
    flt m=ldexp(l+r,-1);
    flt fl=f(l),fm=f(m),fr=f(r);
    return adaptive(l,r,m,eps*15,(r-l)/6*(fl+fr+ldexp(fm,2)),fl,fr,fm,maxrec);
}
}}}

typedef double flt;
flt adaptive(flt l,flt r,flt m,flt eps,flt s,flt fl,flt fr,flt fm,int lim){
    flt h=(r-l)/12;
    flt u=ldexp(l+m,-1),v=ldexp(m+r,-1);
    flt fu=f(u),fv=f(v);
    flt sl=(fl+fm+ldexp(fu,2))*h,
        sr=(fm+fr+ldexp(fv,2))*h,
        s2=sl+sr;
    if(lim<=0 || fabs(s2-s)<= eps)
        return s2+(s2-s)/15;
    eps=ldexp(eps,-1);--lim;
    return  adaptive(l,m,u,eps,sl,fl,fm,fu,lim)
        +   adaptive(m,r,v,eps,sr,fm,fr,fv,lim);
}
flt simpson(flt l,flt r,flt eps=1e-6,int maxrec=10){
    flt m=ldexp(l+r,-1);
    flt fl=f(l),fm=f(m),fr=f(r);
    return adaptive(l,r,m,eps*15,(r-l)/6*(fl+fr+ldexp(fm,2)),fl,fr,fm,maxrec);
}