CF618G Combining Slimes
输出是浮点数,发现合成到 \(50\) 以上的数字的概率已经很小了,可以忽略。
设 \(a_{i,j},b_{i,j}\) 表示用长度为 \(i\) 的格子合成数字 \(j\) 的概率,其中 \(b_{i,j}\) 表示第一个数字必须为 \(2\),得 \(a_{i,j}=a_{i,j-1}\times a_{i-1,j-1},b_{i,j}=b_{i,j-1}\times a_{i-1,j-1}\)。不难发现当 \(i>j\) 时有 \(a_{i+1,j}=a_{i,j},b_{i+1,j}=b_{i,j}\)。
设 \({a}'_{i,j},{b}'_{i,j}\) 表示最终序列倒数第 \(i\) 个格子合成数字 \(j\) 的概率,得 \({a}'_{i,j}=a_{i,j}(1-a_{i-1,j}),{b}'_{i,j}=b_{i,j}(1-a_{i-1,j})\)。
考虑 \(DP\),设 \(f_{i,j}\) 表示当最终序列倒数第 \(i\) 个格子为 \(j\) 时,倒数 \(i\) 个数字和的期望,得:
\[\large f_{i,j}=
\begin{cases}
\frac{1}{\sum\limits_{k=1}^{j-1}{a}'_{i-1,k}}\left(\sum\limits_{k=1}^{j-1}f_{i-1,k}\times {a}'_{i-1,k} \right)+j &j\neq 1 \\
\frac{1}{\sum\limits_{k=2}^{50}{b}'_{i-1,k}}\left(\sum\limits_{k=2}^{50}f_{i-1,k}\times {b}'_{i-1,k} \right)+j &j=1
\end{cases}
\]
预处理前 \(50\) 项后矩阵快速幂即可。
#include<bits/stdc++.h>
#define maxn 60
using namespace std;
template<typename T> inline void read(T &x)
{
x=0;char c=getchar();bool flag=false;
while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
if(flag)x=-x;
}
int n,lim=50;
double p,ans,s;
double a[maxn][maxn],b[maxn][maxn],f[maxn][maxn];
struct matrix
{
double a[maxn][maxn];
}m,v;
matrix operator * (const matrix &x,const matrix &y)
{
matrix z;
memset(z.a,0,sizeof(z.a));
for(int k=0;k<=lim;++k)
for(int i=0;i<=lim;++i)
for(int j=0;j<=lim;++j)
z.a[i][j]+=x.a[i][k]*y.a[k][j];
return z;
}
int main()
{
read(n),scanf("%lf",&p),p/=(double)1000000000;
for(int i=1;i<=lim;++i) a[i][1]=p,a[i][2]=b[i][2]=1-p;
for(int i=2;i<=lim;++i)
for(int j=2;j<=lim;++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=lim;i;--i)
for(int j=1;j<=lim;++j)
a[i][j]*=1-a[i-1][j],b[i][j]*=1-a[i-1][j];
f[1][1]=1,f[1][2]=2;
for(int i=2;i<=lim;++i)
{
for(int j=2;j<=lim;++j)
{
s=0;
for(int k=1;k<j;++k) f[i][j]+=f[i-1][k]*a[i-1][k],s+=a[i-1][k];
f[i][j]=f[i][j]/s+j;
}
s=0;
for(int k=2;k<=lim;++k) f[i][1]+=f[i-1][k]*b[i-1][k],s+=b[i-1][k];
f[i][1]=f[i][1]/s+1;
}
if(n<=lim)
{
for(int i=1;i<=n+1;++i) ans+=f[n][i]*a[n][i];
printf("%.15lf",ans);
return 0;
}
v.a[0][0]=m.a[0][0]=1;
for(int i=1;i<=lim;++i) v.a[0][i]=f[lim][i],m.a[0][i]=i;
for(int i=2;i<=lim;++i)
{
s=0;
for(int j=1;j<i;++j) m.a[j][i]+=a[lim][j],s+=a[lim][j];
for(int j=1;j<i;++j) m.a[j][i]/=s;
}
s=0;
for(int j=2;j<=lim;++j) m.a[j][1]+=b[lim][j],s+=b[lim][j];
for(int j=2;j<=lim;++j) m.a[j][1]/=s;
n-=lim;
while(n)
{
if(n&1) v=v*m;
m=m*m,n>>=1;
}
for(int i=1;i<=lim;++i) ans+=v.a[0][i]*a[lim][i];
printf("%.15lf",ans);
return 0;
}