数值积分

template<class T = double>
struct Simpson {
    function<T(T)> f;
    Simpson(function<T(T)> f) : f(f) {}
    T simpson(T L, T R) {
        T mid = (L + R) / 2;
        return (R - L) * (f(L) + 4 * f(mid) + f(R)) / 6;
    }
    T asr(T L, T R, T eps, T ans) {
        T mid = (L + R) / 2;
        T fL = simpson(L, mid), fR = simpson(mid, R);
        if (fabs(fL + fR - ans) <= 15 * eps) {
            return fL + fR + (fL + fR - ans) / 15;
        }
        return asr(L, mid, eps / 2, fL) + asr(mid, R, eps / 2, fR);
    }
    T operator()(T L, T R) {
        return (*this)(L, R, 1e-12);
    }
    T operator()(T L, T R, T eps) {
        return asr(L, R, eps, simpson(L, R));
    }
}; // Simpson

template<class T = double>
struct Romberg {
    function<T(T)> f;
    Romberg(function<T(T)> f) : f(f) {}
    T ctqf(T L, T R, T h) {
        static T ans;
        ans = 0.0;
        for (T i = L + h / 2 ; i < R ; i += h) {
            ans += f(i);
        }
        return ans;
    }
    T operator()(T L, T R, T eps) {
        static T T1, T2, S1, S2, C1, C2, R1, R2, h, S;
        static int cnt, k;
        k = 1; cnt = 0; h = R - L;
        T1 = (f(L) + f(R)) * h / 2;
        while (true) {
            S = ctqf(L, R, h); cnt++;
            T2 = (T1 + h * S) / 2; S2 = (4 * T2 - T1) / 3; h /= 2;
            T1 = T2; S1 = S2; C1 = C2; R1 = R2;
            if (k == 1) {
                k++;
                continue;
            }
            C2 = (16 * S2 - S1) / 15;
            if (k == 2) {
                k++;
                continue;
            }
            R2 = (64 * C2 - C1) / 63;
            if (k == 3) {
                k++;
                continue;
            }
            if (fabs(R1 - R2) < eps || cnt > 100) {
                break;
            }
        }
        return R2;
    }
}; // Romberg