题目描述:
洛谷 link
1
0
−
9
<
p
≤
1
10^{-9}<p\le 1
10−9<p≤1
题目分析:
先不考虑格子长度的限制,我们要得到数字 x x x 的概率 p [ x ] = p [ x − 1 ] 2 p[x]=p[x-1]^2 p[x]=p[x−1]2,而 p [ 2 ] m a x = 1 − 1 0 − 9 p[2]_{max}=1-10^{-9} p[2]max=1−10−9,那么 p [ x ] m a x = ( 1 − 1 0 − 9 ) 2 x − 1 p[x]_{max}=(1-10^{-9})^{2^{x-1}} p[x]max=(1−10−9)2x−1,当 x = 50 x=50 x=50 时已经无限趋于0,可以近似认为不可能得到大于 50 的数。
设
a
[
i
]
[
j
]
a[i][j]
a[i][j] 表示用
i
i
i 个格子得到第一个数为
j
j
j 的概率,有:
{
j
>
2
,
a
[
i
]
[
j
]
=
a
[
i
]
[
j
−
1
]
∗
a
[
i
−
1
]
[
j
−
1
]
a
[
1
]
[
1
]
=
p
,
a
[
1
]
[
2
]
=
1
−
p
i
>
1
,
a
[
i
]
[
1
]
=
1
,
a
[
i
]
[
2
]
=
1
−
p
+
p
2
\begin{cases} j>2, ~a[i][j]=a[i][j-1]*a[i-1][j-1]\\ a[1][1]=p,~a[1][2]=1-p\\ i>1,~a[i][1]=1,~a[i][2]=1-p+p^2\end{cases}
⎩⎪⎨⎪⎧j>2, a[i][j]=a[i][j−1]∗a[i−1][j−1]a[1][1]=p, a[1][2]=1−pi>1, a[i][1]=1, a[i][2]=1−p+p2
如果相邻两个格子为
1
2
1~2
1 2,那么就不会往左合并了,可能会用到第一个产生的数为 2 形成
j
j
j 的概率。
设
b
[
i
]
[
j
]
b[i][j]
b[i][j] 表示用
i
i
i 个格子,第一次生成的数为2,得到数
j
j
j 的概率,有:
{
j
>
2
,
b
[
i
]
[
j
]
=
b
[
i
]
[
j
−
1
]
∗
a
[
i
−
1
]
[
j
−
1
]
b
[
i
]
[
2
]
=
1
−
p
\begin{cases}j>2,~b[i][j]=b[i][j-1]*a[i-1][j-1]\\ b[i][2]=1-p\end{cases}
{j>2, b[i][j]=b[i][j−1]∗a[i−1][j−1]b[i][2]=1−p
可以发现当 i > j i>j i>j 时, a [ i + 1 ] [ j ] = a [ i ] [ j ] , b [ i + 1 ] [ j ] = b [ i ] [ j ] a[i+1][j]=a[i][j],~b[i+1][j]=b[i][j] a[i+1][j]=a[i][j], b[i+1][j]=b[i][j]
设 d p [ i ] [ j ] dp[i][j] dp[i][j] 表示在最终序列中从后往前数第 i i i 个数为 j j j 的前提下,后 i i i 个数的期望大小和。
对于
j
>
1
j>1
j>1,转移时要枚举第
i
−
1
i-1
i−1 个数为
k
<
j
k<j
k<j,考虑在最终序列中第
i
i
i 个数为
j
j
j 的前提下,第
i
−
1
i-1
i−1 个数为
k
k
k 的概率,不难发现它等于
a
[
i
−
1
]
[
k
]
∗
(
1
−
a
[
i
−
1
]
[
k
]
)
∑
s
=
1
j
−
1
a
[
i
−
1
]
[
s
]
∗
(
1
−
a
[
i
−
1
]
[
s
]
)
\frac {a[i-1][k]*(1-a[i-1][k])}{\sum_{s=1}^{j-1}a[i-1][s]*(1-a[i-1][s])}
∑s=1j−1a[i−1][s]∗(1−a[i−1][s])a[i−1][k]∗(1−a[i−1][k])
因为是在最终序列中,所以必须要保证“稳定”,因此记
A
[
i
]
[
j
]
=
a
[
i
]
[
j
]
∗
(
1
−
a
[
i
−
1
]
[
j
]
)
A[i][j]=a[i][j]*(1-a[i-1][j])
A[i][j]=a[i][j]∗(1−a[i−1][j]) 表示第
i
i
i 个数为
j
j
j 且第
i
−
1
i-1
i−1 个数
<
j
<j
<j 的概率。(特别的,
A
[
2
]
[
1
]
A[2][1]
A[2][1] 的第二个数是可以为 2 的,因此不等于0)
所以对于
j
>
1
j>1
j>1,有:
d
p
[
i
]
[
j
]
=
j
+
∑
k
=
1
j
−
1
d
p
[
i
−
1
]
[
k
]
∗
A
[
i
−
1
]
[
k
]
∑
s
=
1
j
−
1
A
[
i
−
1
]
[
s
]
dp[i][j]=j+\sum_{k=1}^{j-1}dp[i-1][k]*\frac {A[i-1][k]}{\sum_{s=1}^{j-1}A[i-1][s]}
dp[i][j]=j+k=1∑j−1dp[i−1][k]∗∑s=1j−1A[i−1][s]A[i−1][k]
而当
j
=
1
j=1
j=1 时,要枚举第
i
−
1
i-1
i−1 个数为
k
>
1
k>1
k>1,第一个生成的数必须是2,同理为
k
k
k 的概率为
b
[
i
−
1
]
[
k
]
∗
(
1
−
a
[
i
−
1
]
[
k
]
)
∑
s
=
2
50
b
[
i
−
1
]
[
s
]
∗
(
a
[
i
−
1
]
[
s
]
)
\frac {b[i-1][k]*(1-a[i-1][k])}{\sum_{s=2}^{50}b[i-1][s]*(a[i-1][s])}
∑s=250b[i−1][s]∗(a[i−1][s])b[i−1][k]∗(1−a[i−1][k]),记
B
[
i
]
[
j
]
=
b
[
i
]
[
j
]
∗
(
1
−
a
[
i
−
1
]
[
j
]
)
B[i][j]=b[i][j]*(1-a[i-1][j])
B[i][j]=b[i][j]∗(1−a[i−1][j]),则有:
d
p
[
i
]
[
1
]
=
1
+
∑
k
=
2
50
d
p
[
i
−
1
]
[
k
]
∗
B
[
i
−
1
]
[
k
]
∑
s
=
2
50
B
[
i
−
1
]
[
s
]
dp[i][1]=1+\sum_{k=2}^{50}dp[i-1][k]*\frac {B[i-1][k]}{\sum_{s=2}^{50}B[i-1][s]}
dp[i][1]=1+k=2∑50dp[i−1][k]∗∑s=250B[i−1][s]B[i−1][k]
当 i > j i>j i>j 时 A , B A,B A,B 不会随 i i i 改变,所以先暴力算出 d p [ 51 ] [ . . . ] dp[51][...] dp[51][...] (实际上dp[50]也没问题)后用 51*51 的矩阵加速DP就可以了。
最后 A n s = ∑ j = 1 50 A [ n ] [ j ] ∗ d p [ n ] [ j ] Ans=\sum_{j=1}^{50}A[n][j]*dp[n][j] Ans=∑j=150A[n][j]∗dp[n][j]
Code:
#include<bits/stdc++.h>
#define maxn 55
using namespace std;
const int N = 50;
int n;
double p,a[maxn][maxn],b[maxn][maxn],A[maxn][maxn],B[maxn][maxn],f[maxn][maxn],ans;
struct Mat{
double s[maxn][maxn];
Mat(){memset(s,0,sizeof s);}
double* operator [] (int i){return s[i];}//awesome!
Mat operator * (const Mat &B)const{
Mat ret;
for(int k=0;k<=N;k++) for(int i=0;i<=N;i++) if(s[i][k]) for(int j=0;j<=N;j++)
ret.s[i][j]+=s[i][k]*B.s[k][j];
return ret;
}
}F,G;
int main()
{
scanf("%d%lf",&n,&p),p/=1e9;
a[1][1]=p,a[1][2]=b[1][2]=1-p;
for(int i=2;i<=N;i++){
a[i][1]=p,a[i][2]=1-p+p*p,b[i][2]=1-p;
for(int j=3;j<=i+1;j++) a[i][j]=a[i][j-1]*a[i-1][j-1],b[i][j]=b[i][j-1]*a[i-1][j-1];
}
for(int i=1;i<=N;i++) for(int j=1;j<=i+1;j++) A[i][j]=a[i][j]*(1-a[i-1][j]),B[i][j]=b[i][j]*(1-a[i-1][j]);
for(int i=1;i<=N;i++) f[1][i]=i;
for(int i=2;i<=N;i++){
double s=0;
for(int j=2;j<=N;j++) f[i][1]+=f[i-1][j]*B[i-1][j],s+=B[i-1][j];
f[i][1]=f[i][1]/s+1;
double x=0,y=0;
for(int j=2;j<=N;j++) x+=f[i-1][j-1]*A[i-1][j-1],y+=A[i-1][j-1],f[i][j]=j+x/y;
}
if(n<=N){
for(int j=1;j<=N;j++) ans+=A[n][j]*f[n][j];
return printf("%.8f\n",ans),0;
}
F[0][0]=G[0][0]=1;
for(int i=1;i<=N;i++) G[0][i]=i,F[0][i]=f[N][i];
double s=0;
for(int i=2;i<=N;i++) s+=B[N][i];
for(int k=2;k<=N;k++) G[k][1]=B[N][k]/s;
s=0;
for(int j=2;j<=N;j++){
s+=A[N][j-1];
for(int k=1;k<j;k++) G[k][j]=A[N][k]/s;
}
n-=N;
for(;n;n>>=1,G=G*G) if(n&1) F=F*G;
for(int j=1;j<=N;j++) ans+=A[N][j]*F[0][j];
printf("%.8f\n",ans);
}