SDOI2017 龙与地下城

$$ \text{题意:} \\ \text{若干个随机变量}x_1,x_2…x_Y \\ \text{每个随机变量在}\left[ \text{0,}X-1 \right] \text{内等概率随机取值} \\ \text{求}\sum_{1\le i\le Y}{x_i}\text{落在区间}\left[ A,B \right] \text{内的概率} \\ \text{考虑概率生成函数}F\left( x \right) =\sum_{0\le i\le X-1}{x^i\times \frac{1}{X}} \\ \left( F\left( x \right) \right) ^Y\text{中,}k\text{次项前的系数即为这些随机变量的和为}k\text{的概率。} \\ \,\,\text{把次数在}\left[ A,B \right] \text{之间的系数加起来即可} \\ \text{但这样的复杂度是}O\left( XY\log n \right) \text{的} \\ \text{考虑中心极限定理} \\ \text{设随机变量}\overline{x}=\frac{\sum_{1\le i\le Y}{x_i}}{Y} \\ \text{每个随机变量}x_i\text{的期望为}a=\frac{X-1}{2},\text{方差为}\sigma ^2=\frac{X^2-1}{12} \\ \text{在}Y\text{足够大的时候，可以近似认为}\overline{x}\text{服从正态分布}N\left( a,\frac{\sigma ^2}{Y} \right) \\ \text{所以大力辛普森积分}\int_{\frac{A}{Y}}^{\frac{B}{Y}}{\frac{e^{-\frac{\left( x-a \right) ^2}{2\sigma ^2}}}{\sqrt{2\pi \sigma ^2}}}dx,\text{其中}a=\frac{X-1}{2},\sigma ^2=\frac{X^2-1}{12Y} \\ \text{但是注意对于一些正态分布函数,其在大部分位置取值为0，只在少部分位置取值非0,这时候我们} \\ \text{可以强行让辛普森积分函数递归若干层来保证结果正确} \\ $$

#include <assert.h>
#include <algorithm>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
using std::cin;
using std::cout;
using std::endl;
using int_t = long long int;
using real_t = long double;
using cpx_t = std::complex<real_t>;
const real_t PI = std::acos(-1);
const int_t LARGE = 1 << 22;
const real_t EPS = 1e-5;
template <int_t arg = 1>
void transform(cpx_t* A, int_t size);
int_t upper2n(int_t x) {
    int_t res = 1;
    while (res < x) res *= 2;
    return res;
}
template <class Func>
real_t S(real_t a, real_t b, Func f) {
    return (b - a) / 6 * (f(a) + 4 * f((a + b) / 2) + f(b));
}
template <class Func>
real_t Simpson(real_t a, real_t b, Func f, int_t depth = 0) {
    real_t mid = (a + b) / 2;
    if (std::abs(S(a, b, f) - S(a, mid, f) - S(mid, b, f)) / 15 < EPS &&
        depth > 10)
        return S(a, b, f);
    return Simpson(a, mid, f, depth + 1) + Simpson(mid, b, f, depth + 1);
}
void process1(real_t x, real_t y) {
    static cpx_t A[LARGE + 1];
    static real_t prefix[LARGE + 1];
    int_t size = upper2n(x * y);
    std::fill(A, A + size, cpx_t(0));
    for (int_t i = 0; i <= x - 1; i++) A[i] = (real_t)1 / x;
    transform(A, size);
    for (int_t i = 0; i < size; i++) {
        A[i] = std::pow(A[i], y);
    }
    transform<-1>(A, size);
    for (int_t i = 0; i < size; i++) A[i] /= size;
    prefix[0] = A[0].real();
    for (int_t i = 1; i <= x * y; i++) prefix[i] = prefix[i - 1] + A[i].real();
    for (int_t i = 1; i <= 10; i++) {
        int_t a, b;
        cin >> a >> b;
        if (a > 0)
            cout << prefix[b] - prefix[a - 1] << endl;
        else
            cout << prefix[b] << endl;
    }
}
void process2(real_t x, real_t y) {
    real_t a = (x - 1) / 2;
    real_t sigmap2 = (x * x - 1) / 12 / y;
    for (int_t i = 1; i <= 10; i++) {
        real_t A, B;
        cin >> A >> B;
        cout << Simpson(A / y + EPS, B / y + EPS,
                        [=](real_t x) -> real_t {
                            return expl(-(x - a) * (x - a) / (2 * sigmap2)) /
                                   sqrtl(2 * M_PI * sigmap2);
                        })

             << endl;
    }
}
void process() {
    int_t x, y;
    cin >> x >> y;
    if (x * y <= 3e5)
        process1(x, y);
    else
        process2(x, y);
}
int main() {
    int_t T;
    cin >> T;
    for (int_t i = 1; i <= T; i++) {
        cout.setf(std::ios::fixed);
        cout << std::setprecision(10);
        process();
    }
    return 0;
}
int_t bitRev(int_t x, int_t size2) {
    int_t result = 0;
    for (int_t i = 1; i < size2; i++) {
        result |= (x & 1);
        x >>= 1;
        result <<= 1;
    }
    return result | (x & 1);
}
template <int_t arg = 1>
void transform(cpx_t* A, int_t size) {
    int_t size2 = log2(size);
    for (int_t i = 0; i < size; i++) {
        int_t x = bitRev(i, size2);
        if (x > i) std::swap(A[i], A[x]);
    }
    for (int_t i = 2; i <= size; i *= 2) {
        auto mr = cpx_t(cosl(2 * PI / i), arg * sinl(2 * PI / i));
        for (int_t j = 0; j < size; j += i) {
            cpx_t curr(1, 0);
            for (int_t k = 0; k < i / 2; k++) {
                auto u = A[j + k];
                auto t = curr * A[j + k + i / 2];
                A[j + k] = u + t;
                A[j + k + i / 2] = u - t;
                curr *= mr;
            }
        }
    }
}