[LOJ2267][SDOI2017]龍與地下城-FFT-自適應辛普森積分
龍與地下城
Description
小Q同學是一個熱愛學習的人，但是他最近沉迷於各種遊戲，龍與地下城就是其中之一。
在這個遊戲中，很多場合需要通過擲骰子來產生隨機數，並由此決定角色未來的命運，因此骰子堪稱該遊戲的標誌性道具。
骰子也分為許多種類，比如4面骰、6面骰、8面骰、12面骰、20面骰，其中20面骰用到的機會非常多。當然，現在科技發達，可以用一個隨機數生成器來取代真實的骰子，所以這裡認為骰子就是一個隨機數生成器。
在戰鬥中，骰子主要用來決定角色的攻擊是否命中，以及命中後造成的傷害值。舉個例子，假設現在已經確定能夠命中敵人，那麼
眾所周知，骰子顯示每個數的概率應該是相等的，也就是說，對於一個XX面骰子，顯示
更形式地說，這個骰子顯示的數W滿足離散的均勻分佈，其分佈列為
除此之外還有一些性質
W的一階原點矩(期望)為
W的二階中心矩(方差)為
言歸正傳，現在小Q同學面對著一個生命值為A的沒有防禦的敵人，能夠發動一次必中的
因為小Q同學非常謹慎，他會進行10次模擬戰，每次給出敵人的生命值A以及overkill的標準B，他想知道此時取得屬於他的勝利的概率是多少，你能幫幫他嗎？
Input
第一行是一個正整數T，表示測試資料的組數，
對於每組測試資料：
第一行是兩個整數X , Y，分別表示骰子的面數以及骰子的個數；
接下來10行，每行包含兩個整數A , B，分別表示敵人的生命值A以及overkill的標準B。
Output
對於每組測試資料，輸出10行，對每個詢問輸出一個實數，要求絕對誤差不超過
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
對於
被題目名字吸引進來了……
居然沒有找到題解……
然後就被坑了好久
思路:
很顯然可以有一個DP。
令
那麼顯然很好轉移:
然後可以發現這是個卷積，那麼FFT成點值表示式後快速冪直接上即可~
才怪。這隻有70分。
正當咱完全沒有思路的時候，知乎上此題的出題人說，這題的idea 基於“中心極限定理”。
學習一下這個東西可以發現，其實出題人想要讓咱們用正態分佈：
所謂正態分佈，就是滿足位置引數為
其中，
具體詳見百度百科，本蒟蒻水平有限~
可以發現出題人已經在題面中將期望和方差給出了。然而並沒有什麼人知道這兩句話有用
同樣可以發現然而咱還是沒發現，這題的“擲骰子”事件滿足正態分佈。
那麼咱們可以根據正態分佈的曲線來估計答案!
不加證明的（不會證明）給出正態分佈的解析式:
同時給出正態分佈的影象:
然後現在目標轉換為，求圖上兩點間函式影象的面積。
考慮到求函式的面積，還是這麼奇怪的面積，顯然需要估算了~
聽說考場上有人用拆成小梯形的方法過了，然而明顯這樣做精度比較令人著急……
那麼考慮更靠譜一點的方法:自適應辛普森積分。
辛普森積分嘛，就是用一根迷之二次函式曲線，試圖湊出某個奇怪的函式在某一小部分內圍成的面積的方法~原理自行百度~
自適應辛普森積分就是，不停二分割槽間，如果到某一時刻，用辛普森積分對當前的
這演算法的原理有種自欺欺人的感覺
然而這樣積分出面積，只有80pts。
因為剛才咱們使用的是普通的正態分佈，有很多除法，不夠精確。
需要化成標準正態分佈以提高精度。
具體做法是，對於傳入的引數
把
同時將正態分佈的公式換成標準正態分佈:
這樣就少了很多除法，可過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);