输出是浮点数，发现合成到 (50) 以上的数字的概率已经很小了，可以忽略。
设 (a_{i,j},b_{i,j}) 表示用长度为 (i) 的格子合成数字 (j) 的概率，其中 (b_{i,j}) 表示第一个数字必须为 (2)，得 (a_{i,j}=a_{i,j-1} imes a_{i-1,j-1},b_{i,j}=b_{i,j-1} imes 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}=
egin{cases}
frac{1}{sumlimits_{k=1}^{j-1}{a}'_{i-1,k}}left(sumlimits_{k=1}^{j-1}f_{i-1,k} imes {a}'_{i-1,k}
ight)+j &j
eq 1 \
frac{1}{sumlimits_{k=2}^{50}{b}'_{i-1,k}}left(sumlimits_{k=2}^{50}f_{i-1,k} imes {b}'_{i-1,k}
ight)+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;
}