洛谷P3779 [SDOI2017]龙与地下城（概率论+Simpson+FFT）

题面

传送门

题解

orz shadowice

正态分布

正态分布是随机变量(X)的一种概率分布形式。它用一个期望(mu)和方差(sigma^2)就可以描述，记为(N(mu,sigma^2))。
若随机变量(X)服从一个数学期望为(mu)、方差为(sigma^2)的正态分布，记作(X)∼(N(mu,sigma^2))，读作(X)服从(N(mu,sigma^2))。
当(mu=0,sigma=1)时的正态分布称为标准正态分布。

概率密度函数

用来描述连续型随机变量的分布情况。随机变量的取值落在某个区域内的概率，为概率密度函数在该区域的积分。（或者就是(f(x))在该区域内与(x)轴围成的图形面积）
若随机变量(X)∼(N(mu,sigma^2))，则其概率密度函数为

[f(x)=frac{1}{sqrt{2pi}sigma}e^{-frac{(x-mu)^{2}}{2sigma^{2}}} ]

只要证明了一个变量服从正态分布，就可以直接对概率密度函数的这一区间进行积分了

中心极限定理

当样本量n逐渐趋于无穷大时，n个抽样样本的均值的频数逐渐趋于正态分布（无论总体是什么分布）。
该定理说明，设随机变量(X_1,X_2,…,X_n)独立同分布，它们的期望为(mu)、方差为(sigma^2)，当(n)足够大时（满足精确度需求时），随机变量

[Y_{n}=frac{sum_{i=1}^{n}x_{i}-nmu}{sqrt{n}sigma} ]

(Y_n)服从正态分布，求出其范围后就可以直接对正态分布的概率密度函数求积分了。

本题

然后关于本题有

[mu=frac{n-1}{2}，sigma^2=frac{sum_{i=1}^n(i-mu)^2}{n}=frac{n^2-1}{12}\sum_{i=1}^nX_iin[A,B]\Y_nin[frac{A-nmu}{sqrt nsigma},frac{B-nmu}{sqrt nsigma}] ]

那么我们对于(Y_n)的值域直接上自适应辛普森法就行了

然而还有一个问题，就是(n)等于(1)之类的情况我们显然不能认为它满足精度限制

那么我们构造出概率的生成函数，(F(z)=sum_{i=0}^{x-1}{1over x}z^i)，只要求出(F^y(z))的第(l)到(r)项系数就行了。而且这里多项式快速幂可以直接把点值给快速幂

听说余奶奶直接用生成函数艹过去了……完全没用到任何概率论的芝士……orz

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=(1<<19)+5;const double Pi=acos(-1.0),K=1.0/sqrt(2*Pi),eps=1e-7;
int read(){
    R int res,f=1;R char ch;
    while((ch=getchar())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getchar())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
struct cp{
    double x,y;
    cp(){}
    cp(R double xx,R double yy):x(xx),y(yy){}
    inline cp operator +(const cp &b)const{return cp(x+b.x,y+b.y);}
    inline cp operator -(const cp &b)const{return cp(x-b.x,y-b.y);}
    inline cp operator *(const cp &b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
    inline cp operator *(const double &b)const{return cp(x*b,y*b);}
}A[N],O[N];
int r[N],lim,l,X,Y,L;double res;
inline cp ksm(R cp x,R int y){
    R cp res(1,0);
    for(;y;y>>=1,x=x*x)if(y&1)res=res*x;
    return res;
}
void FFT(cp *A,int ty){
    fp(i,0,lim-1)if(i<r[i])swap(A[i],A[r[i]]);
    for(R int mid=1;mid<lim;mid<<=1)
        for(R int j=0;j<lim;j+=(mid<<1))
            fp(k,0,mid-1){
                R cp x=A[j+k],y=A[j+k+mid]*O[mid+k];
                A[j+k]=x+y,A[j+k+mid]=x-y;
            }
    if(ty==-1){
        reverse(A+1,A+lim);
        R double k=1.0/lim;fp(i,0,lim-1)A[i]=A[i]*k;
    }
}
inline double F(R double x){return K*exp(-x*x*0.5);}
inline double simpson(R double l,R double r){return (F(l)+F(r)+4*F((l+r)/2))*(r-l)/6;}
double calc(double l,double r,double eps,double res){
    double mid=(l+r)/2,ql=simpson(l,mid),qr=simpson(mid,r);
    if(fabs(ql+qr-res)<=15*eps)return ql+qr+(ql+qr-res)/15;
    return calc(l,mid,eps/2,ql)+calc(mid,r,eps/2,qr);
}
void solve1(){
    fp(i,0,lim-1)A[i]=cp(0,0),r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(R int i=1;i<lim;i<<=1)fp(k,0,i-1)O[i+k]=cp(cos(Pi*k/i),sin(Pi*k/i));
    R double xx=1.0/X;
    fp(i,0,X-1)A[i].x=xx;
    FFT(A,1);
    fp(i,0,lim-1)A[i]=ksm(A[i],Y);
    FFT(A,-1);
    fp(i,1,lim-1)A[i].x+=A[i-1].x;
    for(R int i=1,l,r;i<=10;++i)l=read(),r=read(),printf("%.7lf
",l?A[r].x-A[l-1].x:A[r].x);
}
void solve2(){
    double l,r,mu=1.0*(X-1)/2,sigma=1.0*(X*X-1)/12,a=mu*Y,b=sqrt(sigma*Y);
    fp(i,1,10){
        l=1.0*(read()-a)/b,r=1.0*(read()-a)/b;
        printf("%.7lf
",calc(0,r,eps,simpson(0,r))-calc(0,l,eps,simpson(0,l)));
    }
}
int main(){
//	freopen("testdata.in","r",stdin);
    for(int T=read();T;--T){
        X=read(),Y=read(),L=X*Y,lim=1,l=0;
        while(lim<=L)lim<<=1,++l;
        lim<N?solve1():solve2();
    }
    return 0;
}