自适应辛普森法求积分

Bridge UVALive - 3485

要计算$\int_{a}^{b}f(x)dx$
将之放到二维坐标系中，就相当于求面积。
三点辛普森公式：$\frac{f(a)+4f(min(a,b))+f(b)}{6}$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;
#define ll long long
#define df double

df a;
df F(df x){ return sqrt(1+4.0*a*a*x*x);}

df simpson(df a,df b)
{
    df c=a+(b-a)/2.0;
    return (F(a)+4*F(c)+F(b))*(b-a)/6.0;
}
df asr(df a,df b,df eps,df A)
{
    df c=a+(b-a)/2.0;
    df 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.0,L)+asr(c,b,eps/2.0,R);
}
df asr(df a,df b,df eps){
    return asr(a,b,eps,simpson(a,b));
}
df parabola(df w,df h)
{
    a=4.0*h/(w*w);
    return asr(0,w/2,1e-5)*2;
}
int main()
{
    int T;
    scanf("%d",&T);
    for(int kase=1;kase<=T;kase++)
    {
        int D,H,B,L;
        scanf("%d %d %d %d",&D,&H,&B,&L);
        int n=(B+D-1)/D;
        df D1=(df)B/n;
        df L1=(df)L/n;
        df x=0,y=H;
        while(y-x>1e-5)
        {
            df m=x+(y-x)/2;
            if(parabola(D1,m)<L1) x=m;
            else y=m;
        }
        if(kase>1) printf("\n");
        printf("Case %d:\n%.2f\n",kase,H-x);
    }
    return 0;
}