CF618G Combining Slimes【条件概率DP】

题目描述：

洛谷 link $10^{-9}<p\le 1$

题目分析：

先不考虑格子长度的限制，我们要得到数字 $x$ 的概率 $p[x]=p[x-1{]}^2$，而 $p[2{]}_{max}=1-10^{-9}$，那么 $p[x{]}_{max}=(1-10^{-9}{)}^{2^{x-1}}$，当 $x=50$ 时已经无限趋于0，可以近似认为不可能得到大于 50 的数。

设 $a[i][j]$ 表示用 $i$ 个格子得到第一个数为 $j$ 的概率，有： $\{\begin{array}{l}j>2,a[i][j]=a[i][j-1]\ast 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{array}$

如果相邻两个格子为$12$，那么就不会往左合并了，可能会用到第一个产生的数为 2 形成 $j$ 的概率。 设 $b[i][j]$ 表示用 $i$ 个格子，第一次生成的数为2，得到数 $j$ 的概率，有： $\{\begin{array}{l}j>2,b[i][j]=b[i][j-1]\ast a[i-1][j-1]\\ b[i][2]=1-p\end{array}$

可以发现当 $i>j$ 时，$a[i+1][j]=a[i][j],b[i+1][j]=b[i][j]$

---

设 $dp[i][j]$ 表示在最终序列中从后往前数第 $i$ 个数为 $j$ 的前提下，后 $i$ 个数的期望大小和。

对于$j>1$，转移时要枚举第 $i-1$ 个数为 $k<j$，考虑在最终序列中第 $i$ 个数为 $j$ 的前提下，第 $i-1$ 个数为 $k$ 的概率，不难发现它等于 $\frac{a[i-1][k]\ast (1-a[i-1][k])}{{\sum }_{s=1}^{j-1}a[i-1][s]\ast (1-a[i-1][s])}$ 因为是在最终序列中，所以必须要保证“稳定”，因此记 $A[i][j]=a[i][j]\ast (1-a[i-1][j])$ 表示第 $i$ 个数为 $j$ 且第 $i-1$ 个数 $<j$ 的概率。（特别的，$A[2][1]$ 的第二个数是可以为 2 的，因此不等于0）

所以对于 $j>1$，有： $$dp[i][j]=j+\sum _{k=1}^{j-1}dp[i-1][k]\ast \frac{A[i-1][k]}{{\sum }_{s=1}^{j-1}A[i-1][s]}$$

而当 $j=1$ 时，要枚举第 $i-1$ 个数为 $k>1$，第一个生成的数必须是2，同理为 $k$ 的概率为 $\frac{b[i-1][k]\ast (1-a[i-1][k])}{{\sum }_{s=2}^{50}b[i-1][s]\ast (a[i-1][s])}$，记 $B[i][j]=b[i][j]\ast (1-a[i-1][j])$，则有： $$dp[i][1]=1+\sum _{k=2}^{50}dp[i-1][k]\ast \frac{B[i-1][k]}{{\sum }_{s=2}^{50}B[i-1][s]}$$

当 $i>j$ 时 $A,B$ 不会随 $i$ 改变，所以先暴力算出 $dp[51][...]$ （实际上dp[50]也没问题）后用 51*51 的矩阵加速DP就可以了。

最后 $Ans={\sum }_{j=1}^{50}A[n][j]\ast 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); }