CF 618G

题目大意

一开始有$N$个空位。给定一个概率$P$。每次你有$P$的概率在第$N$个空位加入一个$1$，有$1 - P$的概率加入一个$2$。加入之后，你会将这个数字往前推，假如当前数字遇到了另外一个和他值一样的数字，设为$v$，那么当前空位的数字变为$v + 1$，并且继续往前推。直到到达边界或遇到一个不一样的数字。游戏结束当且仅当所有空位都放了数字。一种结束局面的权值定义为所有数字的和。问你最终结束时局面的期望权值。

数据范围

$N \leq 10^9$
$10^{-9} \leq P \leq 1$，且$P$的精度范围为$10^{-9}$

题解

设$A_{i,j}$表示当前有$i$个格子，能够凑出数字$j$的概率，注意只能能凑出，没有保证最终一定是。那么很显然的，$A_{i,j} = A_{i,j - 1} * A_{i - 1,j - 1}$。

$B_{i,j}$的意义是一样的。只是保证在一开始放的数字是2。那么
$B_{i,j} = B_{i,j - 1} * A_{i - 1,j - 1}$

那么我们再做些更正，即$A_{i,j},B_{i,j}$的意义变为第$i$个恰好为$j$的概率。那么
$A_{i,j} = A'_{i,j} * (1 - A'_{i - 1,j})$，$A'$是上面求的$A$。

而且可以发现的一点是，当$j$很大时，$A_{i,j}$就趋近于0了。因为他大概是以平方级别变化，那么事实上当$j \geq 50$时就等于0了。
还有一点是，当$i$比较大时，$A_{i,j} = A_{i - 1,j}$，分析方法同理。那么事实上，只需要预处理出$1 \leq i \leq 50$的值就好了。

接来下考虑怎么求答案。

设$dp_{i,j}$表示当前有$i$个格子，第一个格子的值为$j$，且不会与其他数字结合的期望权值。我们对$j$进行分类讨论。

1. $j = 1$
   那么只需要第二个格子不是1就好了。可以得到递推式
   
   $$dp_{i,j} = j + \frac{\sum_{k=2}^{50} dp_{i-1,k} * B_{i - 1,k}}{\sum_{k=2}^{50} B_{i - 1,k}}$$

   这就是为什么我们要预处理出$B$，因为我们第一个塞进来的不能是1，否则就合并了。
   而且递推式的分母其实就是出现这种情况的概率。

   2.$j > 1$
   可以发现的是，第二个格子的值必然比$j$要小，不然的话肯定会发生合并。那么我们就可以得到另一条式子了
   

   $$dp_{i,j} = j + \frac{\sum_{k=1}^{j-1} dp_{i-1,k} * A_{i - 1,k}}{\sum_{k=1}^{j-1} A_{i - 1,k}}$$

   但这题到这里还没结束，因为$N \leq 10^9$!
   但之前说了，由于$j \leq 50$而且当$N > 50$时$A_{N,j} = A_{N - 1,j}$，也就是说，我们只需要求出$dp_{1->50}$，对于$N > 50$的，每次转移的系数都是一样的！那么我们就可以构出转移的矩阵，直接转移就好了。

   那么我们最终的答案就是
   

   $$\sum_{j=1}^{50} dp_{N,j} * A_{N,j}$$

   代码：

   #include <cstdio>
   #include <cstring>
   #include <algorithm>
   
   using namespace std;
   
   const int MAXN = 55;
   
   typedef double Matrix[MAXN][MAXN];
   
   Matrix A,B,Trans,Dp,Cur;
   double Out[55];
   int N;
   
   void Mul(Matrix &a,Matrix &b,Matrix &c)
   {
       static Matrix Tmp;
       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];
       memcpy(c,Tmp,sizeof c);
   }
   
   int main()
   {
   //  freopen("data.in","r",stdin),freopen("data.out","w",stdout);
       double p,q;
       scanf("%d%lf", &N, &p);
       p = p / (1e9),q = 1 - p;
       for(int i = 1;i <= 50;i ++)
           for(int j = 1;j <= 50;j ++)
           {
               if (j == 1) A[i][j] += p;
               if (j == 2) A[i][j] += q,B[i][j] += q;
               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 = 50;i;i --)
           for(int j = 1;j <= 50;j ++)
               A[i][j] *= (1 - A[i - 1][j]),B[i][j] *= (1 - A[i - 1][j]);
       for(int j = 1;j <= 50;j ++) Dp[1][j] = j;
       for(int i = 2;i <= 50;i ++)
           for(int j = 1;j <= 50;j ++)
           {
               double s = 0;
               if (j == 1)
               {
                   for(int k = 2;k <= 50;k ++) Dp[i][j] += Dp[i - 1][k] * B[i - 1][k],s += B[i - 1][k];
               } else
                   for(int k = 1;k < j;k ++) Dp[i][j] += Dp[i - 1][k] * A[i - 1][k],s += A[i - 1][k];
               Dp[i][j] /= s;
               Dp[i][j] += j;
           }
       Trans[0][0] = 1;
       for(int j = 1;j <= 50;j ++) Trans[0][j] = j;
       //For j = 1
       double s = 0;
       for(int k = 2;k <= 50;k ++) s += B[50][k];
       for(int k = 2;k <= 50;k ++) Trans[k][1] += B[50][k] / s;
       //For j not equal 1
       for(int j = 2;j <= 50;j ++)
       {
           s = 0;
           for(int k = 1;k < j;k ++) s += A[50][k];
           for(int k = 1;k < j;k ++) Trans[k][j] += A[50][k] / s;
       }   
       int bk = N;
       if (N > 50)
       {
           memcpy(Cur[0],Dp[50],sizeof Cur);
           Cur[0][0] = 1;
           for(N -= 50;N;N >>= 1)
           {
               if (N & 1) Mul(Cur,Trans,Cur);
               Mul(Trans,Trans,Trans);
           }
           memcpy(Out,Cur[0],sizeof Out);
       } else memcpy(Out,Dp[N],sizeof Out);
       double ans = 0;
       for(int j = 1;j <= 50;j ++) ans += Out[j] * A[bk > 50 ? 50 : bk][j];
       printf("%.15f\n", ans);
       return 0;
   }