ankikong

double work(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 mid=(L+R)/2;
    double x1=work(L,mid),x2=work(mid,R),x3=work(L,R);
    if (std::fabs(x1+x2-x3)<eps) return x1+x2;
    else return simpson(L,mid)+simpson(mid,R);
}