【CF618G】Combining Slimes
题意：一个长度为$1 imes n$的网格，每次从最右侧往里推入一个数字1或2（数字会一直跑到最左边的空格子里），加入1的概率为p，2的概率为1-p。如果新加入的数与其左边的那个数相同，都=x，则将二者合并变成x+1。然后继续判断是否能与左边合并（跟2048差不多）。问你当最后格子满时，整个网格中所有数的和的期望值。
$nle 10^9$
题解：cf怎么总喜欢利用浮点数精度来出题啊？！（现在遍地都是模意义下的期望mod 998244353，这个性质完全不敢用）
我们设a[i][j]表示用一个长度为i的网格造出一个数字j的概率，显然有$a[i][j]=a[i][j-1] imes a[i-1][j-1]$（j=1,2特殊处理）。以及b[i][j]表示在a的基础之上，要求第一个数字为2的概率，显然$b[i][j]=b[i][j-1] imes a[i-1][j-1]$。
观察一番你会发现，当j增大时a[i][j]急剧减小，当j=50时便可以忽略不计，所以我们就可以只考虑出现1..50的概率了。并且a[..][j]自从不为0开始便保持不变，于是用a[50][..]完全可以代表a[..][..]。
那么如何表示右数第i个格子的数正好为j的概率呢？不难发现它等于$a[i][j] imes(1-a[i-1][j])$。于是令$a'[i][j]=a[i][j] imes(1-a[i-1][j]),b'[i][j]=b[i][j] imes(1-a[i-1][j])$。
然后就可以DP了，用f[i][j]表示第i个格子是j时，右面i个格子总和的期望值。容易得到DP方程：
$j eq1：f[i][j]=j+{(sumlimits_{k=1}^{j-1}f[i-1][k] imes a'[i-1][k])over sumlimits_{k=1}^{j-1}a'[i-1][k]}$
$j=1：f[i][j]=j+{(sumlimits_{k=2}^{50}f[i-1][k] imes b'[i-1][k])over sumlimits_{k=2}^{50}b'[i-1][k]}$
所以先预处理出前50项，剩余的进行矩阵乘法就好了。
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef double db;
const int m=50;
struct M
{
db v[60][60];
M () {memset(v,0,sizeof(v));}
db * operator [] (const int &a) {return v[a];}
M operator * (const M &a) const
{
M b;
int i,j,k;
for(i=0;i<=m;i++) for(j=0;j<=m;j++) for(k=0;k<=m;k++) b.v[i][j]+=v[i][k]*a.v[k][j];
return b;
}
}S,T;
db p,a[60][60],b[60][60],f[60][60];
int n;
inline void pm(int y)
{
while(y)
{
if(y&1) S=S*T;
T=T*T,y>>=1;
}
}
int main()
{
int i,j,k,t;
double tmp;
scanf("%d%d",&n,&t),p=t*1e-9;
a[1][1]=p,a[1][2]=1-p;
b[1][2]=1-p;
for(i=2;i<=m;i++)
{
a[i][1]=p,a[i][2]=1-p,b[i][2]=1-p;
for(j=2;j<=m;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(i=m;i>=1;i--) for(j=1;j<=m;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(i=2;i<=m;i++)
{
for(j=2;j<=m;j++)
{
tmp=0;
for(k=1;k<j;k++) f[i][j]+=f[i-1][k]*a[i-1][k],tmp+=a[i-1][k];
f[i][j]=f[i][j]/tmp+j;
}
tmp=0;
for(k=2;k<=m;k++) f[i][1]+=f[i-1][k]*b[i-1][k],tmp+=b[i-1][k];
f[i][1]=f[i][1]/tmp+1;
}
if(n<=m)
{
tmp=0;
for(i=1;i<=n+1;i++) tmp+=f[n][i]*a[n][i];
printf("%.12lf",tmp);
return 0;
}
S[0][0]=T[0][0]=1;
for(i=2;i<=m;i++)
{
tmp=0;
for(j=1;j<i;j++) T[j][i]+=a[m][j],tmp+=a[m][j];
for(j=1;j<i;j++) T[j][i]/=tmp;
T[0][i]=i;
}
tmp=0;
for(i=2;i<=m;i++) T[i][1]+=b[m][i],tmp+=b[m][i];
for(i=2;i<=m;i++) T[i][1]/=tmp;
T[0][1]=1;
for(i=1;i<=m;i++) S[0][i]=f[m][i];
pm(n-m);
tmp=0;
for(i=1;i<=m;i++) tmp+=S[0][i]*a[m][i];
printf("%.12lf",tmp);
return 0;
}