这题真的太神辣,%了一发题解,原来还能这么搞QWQ
设\(A_{i,j}\)表示不加任何限制时,第\(i\)个格子会出现权值为\(j\)的史莱姆的概率,则有:
\[A_{i,j}=A_{i,j-1}*A_{i-1,j-1}\]
再设\(B_{i,j}\)表示第一个放置的史莱姆权值为\(2\)时(之后放置的史莱姆不管),第\(i\)个格子会出现权值为\(j\)的史莱姆的概率,则有:
\[B_{i,j}=B_{i,j-1}*A_{i-1,j-1}\]
好的现在就可以求第\(i\)个格子刚好会出现权值为\(j\)的史莱姆的概率,分别设它们为\(A'_{i,j}\)和\(B'_{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\)时,右数第\(1\)~\(i\)个格子上的史莱姆权值和,根据期望\(Dp\)的基本套路,我们有:
\[ F_{i,j}=\left\{ \begin{aligned} \frac{ \sum_{k=1}^{j-1} F_{i-1,k}*A'_{i-1,k} }{ \sum_{k=1}^{j-1}A'_{i-1,k} }+j&&(j>1) \\ \\\frac{ \sum_{k=2}^{i+1} F_{i-1,k}*B'_{i-1,k} }{ \sum_{k=2}^{i+1}B'_{i-1,k} }+1&&(j=1) \end{aligned} \right.\]
最终答案为:
\[Ans= \sum\limits_{j=1}^{n+1} A_{n,j}*F_{n,j} \]
于是 丧心病狂的 出题人把\(n\)出到了\(1e9\),但是又很良心地保留12位小数,似乎可以利用一下浮点数精度?
观察\(A\)、\(B\)的递推式,发现它是以平方级别增长的。引用下官方题解:当\(j=50\)时, \(A_{i,j}\approx 10^{-300}\).
所以当\(j\geq 50\)时,可以认为\(A_{i,j}=A_{i,j-1}\),同理,当\(i\geq 50\)时,也有\(A_{i,j}=A_{i-1,j}\)
所以就处理出来所有\(1\leq i,j \leq n\)的\(A'_{i,j},B'_{i,j},F_{i,j}\),剩下的\(F_{i,j}\)转移系数是一样的,直接矩乘优化即可
代码:
#include <bits/stdc++.h>
#define lb double
using namespace std;
lb A[55][55],B[55][55],_A[55][55],_B[55][55],F[55][55];
lb E[55][55],Mar[55][55],tmp[55][55];
void init(int k){
int i,j;
A[1][1]=(lb)k*1e-9,B[1][2]=(1e9-k)*1e-9,A[1][2]=B[1][2];
for(i=2;i<=50;++i){
A[i][1]=A[1][1],B[i][2]=B[1][2],A[i][2]=B[i][2]+A[1][1]*A[1][1];
for(j=3;j<=50;++j){
A[i][j]=A[i-1][j-1]*A[i][j-1];
B[i][j]=B[i][j-1]*A[i-1][j-1];
}
}
_A[1][1]=A[1][1],_A[1][2]=A[1][2],_B[1][2]=B[1][2];
for(i=2;i<=50;++i){
_A[i][1]=A[i][1]*(1.00-A[i-1][1]),_B[i][1]=0;
for(j=2;j<=50;++j){
_A[i][j]=A[i][j]*(1.00-A[i-1][j]);
_B[i][j]=B[i][j]*(1.00-A[i-1][j]);
}
}
for(i=1;i<=50;++i) F[1][i]=i*1.00;
for(i=2;i<=50;++i){
lb x=0,y=0;
for(j=2;j<=50;++j)
x+=F[i-1][j]*_B[i-1][j],y+=_B[i-1][j];
F[i][1]=x/y+1.00;
for(j=2;j<=50;++j){
x=0,y=0;
for(k=1;k<j;++k)
x+=F[i-1][k]*_A[i-1][k],y+=_A[i-1][k];
F[i][j]=x/y+j*1.00;
}
}
E[0][0]=E[0][1]=1;
for(i=2;i<=50;++i){
lb y=0;Mar[0][i]=F[50][i];
for(j=1;j<i;++j) y+=_A[50][j];
for(j=1;j<i;++j) E[j][i]=_A[50][j]/y;
E[0][i]=i;
}
for(i=2;i<=50;++i){
lb y=0;
for(j=2;j<=50;++j) y+=_B[50][j];
E[i][1]=_B[50][i]/y;
}
Mar[0][0]=1;
for(i=1;i<=50;++i) Mar[0][i]=F[50][i];
}
void mul(lb a[][55],lb b[][55]){
memset(tmp,0,sizeof(tmp));
for(int i=0;i<=50;++i)
for(int k=0;k<=50;++k)
for(int j=0;j<=50;++j)
tmp[i][j]+=a[i][k]*b[k][j];
for(int i=0;i<=50;++i)
for(int j=0;j<=50;++j) a[i][j]=tmp[i][j];
}
int main(){
int n,k,i;
lb ans=0;
scanf("%d%d",&n,&k);
init(k);
if(n<=50){
for(i=1;i<=50;++i) ans+=F[n][i]*_A[n][i];
} else{
for(i=n-50;i;i>>=1){
if(i&1) mul(Mar,E);
mul(E,E);
}
for(i=1;i<=50;++i) ans+=Mar[0][i]*_A[50][i];
}
printf("%.12f",ans);
}