自适应Simpson公式

参考刘汝佳<算法指南>P163

#include<cstdio>
#include<cmath>
double a;
double F(double x){
    return sqrt(1+4*a*a*x*x);
}
double simpson(double a,double b)
{
    double c=a+(b-a)/2;
    return (F(a)+4*F(c)+F(b))*(b-a)/6;
}
double asr(double a,double b,double eps,double A)
{
    double c=a+(b-a)/2;
    double L=simpson(a,c),R=simpson(c,b);
    if(fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0;
    return asr(a,c,eps/2,L)+asr(c,b,eps/2,R); 
}
double ASR(double a,double b,double eps){
    return asr(a,b,eps,simpson(a,b));
}
double parabola_arc_length(double w,double h)
{
    a=4.0*h/(w*w);
    return ASR(0,w/2,1e-5)*2;
}
int main()
{
    int T;
    scanf("%d",&T);
    for(int cas=1;cas<=T;cas++){
        int D,H,B,L;
        scanf("%d%d%d%d",&D,&H,&B,&L);
        int n=(B+D-1)/D;
        double D1=(double)B/n;
        double L1=(double)L/n;
        double x=0,y=H;
        while(y-x>1e-5){
            double m=x+(y-x)/2;
            if(parabola_arc_length(D1,m)<L1) x=m;
            else y=m;
        }
        if(cas>1) printf("
");
        printf("Case %d:
%.2lf
",cas,H-x);
    }    
}