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;
}