[LOJ2267][SDOI2017]龙与地下城-FFT-自适应辛普森积分

龙与地下城

Description

小Q同学是一个热爱学习的人，但是他最近沉迷于各种游戏，龙与地下城就是其中之一。

在这个游戏中，很多场合需要通过掷骰子来产生随机数，并由此决定角色未来的命运，因此骰子堪称该游戏的标志性道具。

骰子也分为许多种类，比如4面骰、6面骰、8面骰、12面骰、20面骰，其中20面骰用到的机会非常多。当然，现在科技发达，可以用一个随机数生成器来取代真实的骰子，所以这里认为骰子就是一个随机数生成器。

在战斗中，骰子主要用来决定角色的攻击是否命中，以及命中后造成的伤害值。举个例子，假设现在已经确定能够命中敌人，那么$YdX$（也就是掷出Y个X面骰子之后所有骰子显示的数字之和）就是对敌人的基础伤害。在敌人没有防御的情况下，这个基础伤害就是真实伤害。

众所周知，骰子显示每个数的概率应该是相等的，也就是说，对于一个XX面骰子，显示$0, 1, 2,\dots ,X−1$中每一个数字的概率都是$\frac {1}{x}$

更形式地说，这个骰子显示的数W满足离散的均匀分布，其分布列为

除此之外还有一些性质

W的一阶原点矩(期望)为$v_1(W)=E(W)=\sum_{i=0}^{X-1}iP(W=i)=\frac {X-1}{2}$
W的二阶中心矩(方差)为$\mu_2(W)=E((W-E(W))^2)=\sum_{i=0}^{X-1}(i-E(W))^2P(W=i)=\frac {X^2-1}{12}$

​
言归正传，现在小Q同学面对着一个生命值为A的没有防御的敌人，能够发动一次必中的$YdX$攻击，显然只有造成的伤害不少于敌人的生命值才能打倒敌人。但是另一方面，小Q同学作为强迫症患者，不希望出现overkill，也就是造成的伤害大于B的情况，因此只有在打倒敌人并且不发生overkill的情况下小Q同学才会认为取得了属于他的胜利。

因为小Q同学非常谨慎，他会进行10次模拟战，每次给出敌人的生命值A以及overkill的标准B，他想知道此时取得属于他的胜利的概率是多少，你能帮帮他吗？

Input

第一行是一个正整数T，表示测试数据的组数，
对于每组测试数据：
第一行是两个整数X , Y，分别表示骰子的面数以及骰子的个数；
接下来10行，每行包含两个整数A , B，分别表示敌人的生命值A以及overkill的标准B。

Output

对于每组测试数据，输出10行，对每个询问输出一个实数，要求绝对误差不超过$0.013579$，也就是说，记输出为aa，答案为bb，若满足$|a-b|\leq 0.013579$，则认为输出是正确的。

Sample Input

1
2 19
0 0
0 1
0 2
0 3
0 4
0 5
0 6
0 7
0 8
0 9

Sample Output

0.000002
0.000038
0.000364
0.002213
0.009605
0.031784
0.083534
0.179642
0.323803
0.500000

Hint

对于$100\%$的数据，$T \leq 10$，$2 \leq X \leq 20$，$1 \leq Y \leq 200000$，$0 \leq A \leq B \leq (X-1)Y$，保证满足$Y > 800$的数据不超过22组。

---

被题目名字吸引进来了……
居然没有找到题解……
然后就被坑了好久

---

思路:
很显然可以有一个DP。
令$f[i][j]$表示丢了$i$个骰子，点数为$j$的概率。
那么显然很好转移:
$f[i][j]=\sum_{k \leq j ,k \leq x-1} f[i-1][j-k]*\frac{1}{n}$

然后可以发现这是个卷积，那么FFT成点值表达式后快速幂直接上即可~

才怪。这只有70分。
正当咱完全没有思路的时候，知乎上此题的出题人说，这题的idea 基于“中心极限定理”。

学习一下这个东西可以发现，其实出题人想要让咱们用正态分布：
所谓正态分布，就是满足位置参数为$\mu$、尺度参数为 $\sigma$ 的概率分布的随机变量的分布曲线。
其中，$\mu$为期望，$\sigma$为方差的平方根。
具体详见百度百科，本蒟蒻水平有限~

可以发现出题人已经在题面中将期望和方差给出了。然而并没有什么人知道这两句话有用
同样可以发现然而咱还是没发现，这题的“掷骰子”事件满足正态分布。
那么咱们可以根据正态分布的曲线来估计答案!

不加证明的（不会证明）给出正态分布的解析式:
$f(x)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x-\mu)^2}{2\sigma^2})$
同时给出正态分布的图像:

然后现在目标转换为，求图上两点间函数图像的面积。
考虑到求函数的面积，还是这么奇怪的面积，显然需要估算了~

听说考场上有人用拆成小梯形的方法过了，然而明显这样做精度比较令人着急……
那么考虑更靠谱一点的方法:自适应辛普森积分。

辛普森积分嘛，就是用一根迷之二次函数曲线，试图凑出某个奇怪的函数在某一小部分内围成的面积的方法~原理自行百度~
自适应辛普森积分就是，不停二分区间，如果到某一时刻，用辛普森积分对当前的$(l,r)$积分得到的结果$f(l,r)$，与$f(l,mid)+f(mid,r)$之间的差很小，咱就认为它足够精确了。
这算法的原理有种自欺欺人的感觉

然而这样积分出面积，只有80pts。
因为刚才咱们使用的是普通的正态分布，有很多除法，不够精确。
需要化成标准正态分布以提高精度。

具体做法是，对于传入的参数$x$，进行如下处理:
$y=\frac{x-\mu}{\sigma}$
把$y$作为以前的$x$使用。
同时将正态分布的公式换成标准正态分布:
$f(x)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{x^2}{2})$
这样就少了很多除法，可过100pts!

注意，小数据由于精度误差，仍需要使用dp（咱用的FFT）

#include<bits/stdc++.h>

using namespace std;

inline int read()
{
    int x=0;char ch=getchar();
    while(ch<'0' || '9'<ch)ch=getchar();
    while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();
    return x;
}

typedef double db;

db x,y;

namespace Fast_Fouriet_Transform
{
    typedef complex<db> cp;
    const int M=800009;
    const db pi=acos(-1);

    cp a[M],b[M];
    db per,ans[M];
    int m,l,rev[M];

    inline void FFT(cp *a,int n,int f)
    {
        for(int i=0;i<n;i++)
            if(i<rev[i])
                swap(a[i],a[rev[i]]);
        for(int h=2;h<=n;h<<=1)
        {
            cp w(cos(2*pi/h),f*sin(2*pi/h));
            for(int j=0;j<n;j+=h)
            {
                cp wn(1,0),x,y;
                for(int k=j;k<j+(h>>1);k++)
                {
                    x=a[k],y=wn*a[k+(h>>1)];
                    a[k]=x+y;
                    a[k+(h>>1)]=x-y;
                    wn*=w;
                }
            }
        }
        if(f==-1)
            for(int i=0;i<n;i++)
                a[i].real()/=(db)n;
    }

    int mian()
    {
        per=1.0/(db)x;

        for(m=1,l=0;m<(y*x);m<<=1)l++;
        for(int i=0;i<m;i++)
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));

        memset(a,0,sizeof(a));
        memset(b,0,sizeof(b));

        a[0].real()=1.0;
        for(int i=0;i<x;i++)
            b[i].real()=per;

        FFT(a,m,1);FFT(b,m,1);
        for(int i=0;i<m;i++)
        {
            int cnt=y;
            while(cnt)
            {
                if(cnt&1)
                    a[i]=a[i]*b[i];
                b[i]=b[i]*b[i];
                cnt>>=1;
            }
        }
        FFT(a,m,-1);

        ans[0]=a[0].real();
        for(int i=1;i<m;i++)
            ans[i]=a[i].real()+ans[i-1];

        for(int i=1;i<=10;i++)
        {
            int a=read(),b=read();
            if(a==0)
                printf("%.8lf\n",ans[b]);
            else
                printf("%.8lf\n",ans[b]-ans[a-1]);
        }
        return 0;
    }
}

namespace simpsons
{
    const db pi2=sqrt(acos(-1)*2.0);
    const db eps=1e-12;
    db sigma,mu;
    inline db sqr(db x){return x*x;}
    inline db f(db x){return exp(-sqr(x)/2)/pi2;}
    inline bool dcmp(db x){return fabs(x)<15*eps;}
    inline void init()
    {
        sigma=sqrt((sqr(x)-1.0)/12.0);
        mu=(x-1.0)/2.0;
    }

    inline db simpson(db a,db b,db fl,db fmid,db fr)
    {
        return (b-a)*(fl+4.0*fmid+fr)/6.0;
    }

    inline db solve(db l,db r)
    {
        db mid=(l+r)/2.0,lm=(l+mid)/2.0,mr=(mid+r)/2.0;
        db fl=f(l),fmid=f(mid),fr=f(r),flm=f(lm),fmr=f(mr);

        db sipl=simpson(l,mid,fl,flm,fmid);
        db sipr=simpson(mid,r,fmid,fmr,fr);
        db sipm=simpson(l,r,fl,fmid,fr);

        if(dcmp(sipl+sipr-sipm))
            return sipl+sipr+(sipl+sipr-sipm)/15;
        else return solve(l,mid)+solve(mid,r);
    }

    int mian()
    {
        init();
        db base=solve(0,(x-1)*y);
        for(int i=1;i<=10;i++)
        {
            db a=(((db)read()-0.5)/y-mu)*sqrt(y)/sigma;
            db b=(((db)read()+0.5)/y-mu)*sqrt(y)/sigma;
            printf("%.8lf\n",solve(0,b)-solve(0,a));
        }
    }
}

int main()
{
    int T=read();
    while(T--)
    {
        scanf("%lf%lf",&x,&y);
        if(x*y<=400000)
            Fast_Fouriet_Transform::mian();
        else
            simpsons::mian();
    }

    return 0;
}