Luogu

LOJ

分析

一种 FFT 做法

考虑 DP 。设 $dp_{i,j}$ 表示前 $i$ 个骰子和为 $j$ 的概率。

可以发现转移是一个多项式快速幂，于是直接 FFT 就可以了。

（大概）可以理解成概率生成函数卷积？

然而 FFT 掉精度很严重，我只能过 $60$ 分。所以对于大数据要考虑其它做法。

一种定积分做法

首先要知道中心极限定理的一个推论：

设有 $n$ 个独立同分布的随机变量 $x_1,x_2,\cdots,x_n$ ，满足期望为 $\mu$ ，方差为 $\sigma^2$ 。

则当 $n$ 足够大时，随机变量 $Y_n=\frac{\sum_{i=1}^nx_i-n\mu}{\sqrt{n\sigma^2}}$ 近似地服从正态分布 $N(0,1)$ 。

如果 $\sum_{i=1}^nx_i\in[a,b]$ ，那么 $Y_n\in\left[\frac{a-n\mu}{\sqrt{n\sigma^2}},\frac{b-n\mu}{\sqrt{n\sigma^2}}\right]$ 。

因为 $Y_n$ 近似服从正态分布 $N(0,1)$ ，所以可以知道 $Y_n$ 的概率密度函数为
$$
\varphi(x)=\frac{1}{\sqrt{2\pi}}e^{\frac{-x^2}{2}}
$$
于是直接上自适应辛普森法就行了。注意精度。

代码

// ===================================
//   author: M_sea
//   website: http://m-sea-blog.com/
// ===================================
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define re register
using namespace std;

inline int read() {
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=800000+10;
const double eps=1e-12;

int X,Y;

namespace FFT {
    struct Complex {
        double x,y;
        Complex(double x_=0,double y_=0):x(x_),y(y_) { }
    };
    Complex operator +(Complex a,Complex b) { return Complex(a.x+b.x,a.y+b.y); }
    Complex operator -(Complex a,Complex b) { return Complex(a.x-b.x,a.y-b.y); }
    Complex operator *(Complex a,Complex b) { return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
    inline Complex qpow(Complex a,int b) {
        Complex c(1,0);
        for (;b;b>>=1,a=a*a)
            if (b&1) c=a*c;
        return c;
    }

    int r[N];
    inline void FFT(Complex* a,int n,int op) {
        for (re int i=0;i<n;++i) if (i<r[i]) swap(a[i],a[r[i]]);
        for (re int i=1;i<n;i<<=1) {
            Complex rot(cos(M_PI/i),op*sin(M_PI/i));
            for (re int j=0;j<n;j+=(i<<1)) {
                Complex w(1,0);
                for (re int k=0;k<i;++k,w=w*rot) {
                    Complex x=a[j+k],y=w*a[j+k+i];
                    a[j+k]=x+y,a[j+k+i]=x-y;
                }
            }
        }
        if (op==-1) {
            for (re int i=0;i<n;++i) a[i].x/=n;
        }
    }

    Complex F[N];
    inline void work() {
        int lim=1,l=0;
        for (;lim<=X*Y;lim<<=1,++l);
        for (re int i=1;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        for (re int i=0;i<X;++i) F[i].x=1.0/X;
        for (re int i=X;i<lim;++i) F[i].x=F[i].y=0;
        FFT(F,lim,1);
        for (re int i=0;i<lim;++i) F[i]=qpow(F[i],Y);
        FFT(F,lim,-1);
        for (re int i=1;i<=10;++i) {
            int a=read(),b=read(); double res=0;
            for (re int j=a;j<=b;++j) res+=F[j].x;
            printf("%.6lf\n",res);
        }
    }
}

namespace Simpson {
    inline double f(double x) { return 1/sqrt(2*M_PI)*exp(-x*x/2.0); }
    inline double simpson(double l,double r) {
        double mid=(l+r)/2;
        return (f(l)+4*f(mid)+f(r))*(r-l)/6;
    }
    inline double asr(double l,double r,double ans) {
        double mid=(l+r)/2,L=simpson(l,mid),R=simpson(mid,r);
        if (fabs(L+R-ans)<eps) return ans;
        else return asr(l,mid,L)+asr(mid,r,R);
    }
    inline double calc(double l,double r) { return asr(l,r,simpson(l,r)); }
    inline void work() {
        double mu=(X-1)/2.0,sigma2=(X*X-1)/12.0;
        for (re int i=1;i<=10;++i) {
            double a=read(),b=read();
            a=(a-Y*mu)/sqrt(Y*sigma2),b=(b-Y*mu)/sqrt(Y*sigma2);
            printf("%.6lf\n",calc(0,b)-calc(0,a)); //卡精度
        }
    }
}

int main() {
    int T=read();
    while (T--) {
        X=read(),Y=read();
        if (X*Y<=160000) FFT::work();
        else Simpson::work();
    }
    return 0;
}