自适应辛普森法求积分

Bridge UVALive - 3485

要计算∫baf(x)dx
将之放到二维坐标系中，就相当于求面积。
三点辛普森公式：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 doubledf 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;}