CF618G Combining Slimes 题解

题意：

有一个 $1 \times n$ 的网格，你将对这个网格进行 $n$ 次操作。

每次从右端插入一个数 $1$ 或 $2$ 。当有两个相邻的数字 $x$ 时，他们会合并为一个数字 $x+1$ 。

插入 $1$ 的概率为 $p$ ，插入 $2$ 的概率为 $1-p$ 。（原题给的是整数 $p_0=p\times 10^{9}$ ）

求无法合并时，$n$ 个格子中所有数的和的期望。浮点输出，要求误差小于 $10^{-4}$ 。

输入：一行，$n,p_0$ 。输出：一行，答案。数据范围：$1\le n,p_0 \le 10^9$ 。

由于每个数都可能出现，根据期望的线性性，所以我们可以分别计算每个数的贡献。

设 $a_{i,j}$ 为用 $i$ 个格子且最左边的那一格为 $j$ 的概率，则

$a_{i,j} = a_{i,j-1}\cdot a_{i-1,j-1}$

根据经验，$a_{i,j}$ 可近似为 $(1-p)^{2^{j-1}}$ ，当 $j = 50$ 时这个值最大为 $10^{-10^{5.39}}$ 。

真是小到连高精度浮点数都汗流浃背了$😅$。因此不妨记 $m=50$ ，此后对于 $j > m$ 我们一律视作 $0$ 。

设 $A_{i,j}$ 为用 $i$ 个格子且最左边的那一格为 $j$ 且它不会被合并的概率，则

$A_{i,j} = a_{i,j} \cdot(1 - a_{i - 1, j})$

设 $f_{i,j}$ 为最终序列中右 $i$ 位中最左侧的格子是 $j$ 时，这 $i$ 位的期望和，则

$f_{i,j} = j + \frac{\sum_{k=1}^{j-1}f_{i-1,k}\cdot A_{i-1,k}}{\sum_{k=1}^{j-1}A_{i-1,k}} \quad (j > 1)$

当 $j=1$ 时，下一个插入的可能比 $j$ 大。

不过，注意到下一个插入的一定是 $2$ 。考虑特别处理这种情况。

设 $b_{i,j}$ 为用 $i$ 个格子且最左边的那一格为 $j$ 且第一次插入的是 $2$ 时的概率。

这里转移其实很简单，我们只需要考虑 $2$ 是最先加入的，那么

$b_{i,j} = b_{i,j-1}\cdot a_{i-1,j-1}$

设 $B_{i,j}$ 为用 $i$ 个格子且最左边的那一格为 $j$ 且它不会被合并且第一次插入的是 $2$ 时的概率，则

$B_{i,j} = b_{i,j} \cdot (1 - a_{i-1,j})$

那么

$f_{i,j} = \begin{cases} j + \large\frac{\sum_{k=2}^{m}f_{i-1,k}\, \cdot \,B_{i-1,k}}{\sum_{k=2}^{m}B_{i-1,k}} && j = 1 \\[6pt]j + \large\frac{\sum_{k=1}^{j-1}f_{i-1,k}\, \cdot \,A_{i-1,k}}{\sum_{k=1}^{j-1}A_{i-1,k}} && j > 1 \end{cases}$

最终答案就是

$\sum_{j=1}^m A_{n,j} \cdot f_{n,j}$

由于 $i > m$ 时 $A_{i,j} \approx A_{m,j}$ ，于是答案为

$\sum_{j=1}^m A_{m,j} \cdot f_{n,j}$

那么我们只需要处理出 $f_{n,j}$ 即可。

注意到 $f$ 的转移由上一项递推而来，且 $m = 50$ ，考虑矩阵快速幂

$\small [\begin{array}{l}f_{m,1} & f_{m,2} & f_{m,3} & \cdots &f_{m,m}\end{array}] \times \left[\begin{array}{cccccc} 1 & 1 & 2 & 3 & \cdots & m \\[6pt]0 & 0 & \frac{A_{m, 1}}{\sum_{k=1}^1 A_{m, k}} & \frac{A_{m, 1}}{\sum_{k=1}^2 A_{m, k}} & \cdots & \frac{A_{m, 1}}{\sum_{k=1}^{m-1} A_{m, k}} \\[6pt]0 & \frac{B_{m, 2}}{\sum_{k=2}^m B_{m, k}} & 0 & \frac{A_{m, 2}}{\sum_{k=1}^2 A_{m, k}} & \cdots & \frac{A_{m, 2}}{\sum_{k=1}^{m-1} A_{m, k}} \\[6pt]0 & \frac{B_{m, 3}}{\sum_{k=2}^m B_{m, k}} & 0 & 0 & \cdots & \frac{A_{m, 3}}{\sum_{k=1}^{m-1} A_{m, k}} \\[6pt]\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\[6pt]0 & \frac{B_{m, m}}{\sum_{k=2}^m B_{m, k}} & 0 & 0 & \cdots & 0 \end{array}\right]^{n - m} = [\begin{array}{l}f_{n,1} & f_{n,2} & f_{n,3} & \cdots &f_{n,m}\end{array}]$

时间复杂度 $\mathcal{O}(m^3 \log n)$

代码：

#include <bits/stdc++.h>
using namespace std;
// #define int long long
// #define INF 0x3f3f3f3f3f3f3f3f
typedef long long ll;
typedef long double lf;
typedef vector< vector<lf> > mat;
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define M ((int)(66))

lf a[M][M], b[M][M], c[M][M], d[M][M], dp[M][M];
mat mul(mat &A, mat &B)
{
    int n = A.size(), m = A[0].size(), p = B[0].size();
    mat C(n, vector<lf>(p));
    for(int k = 0; k < m; k++)
        for(int i = 0; i < n; i++)
            for(int j = 0; j < p; j++) C[i][j] += A[i][k] * B[k][j];
    return C;
}
mat qpow(mat A, int k)
{
    int n = max(A.size(), A[0].size());
    mat B(n, vector<lf>(n));
    for(int i = 0; i < n; i++) B[i][i] = 1;
    for(; k; k >>= 1, A = mul(A, A))
        if(k & 1) B = mul(B, A);
    return B;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n; lf p; cin >> n >> p; p /= 1e9; const int m = 50;
    a[1][1] = p; a[1][2] = 1 - p; b[1][2] = 1 - p;
    for(int i = 2; i <= m; i++)
    {
        a[i][1] = p; a[i][2] = 1 - p + p * p; b[i][2] = 1 - p;
        for(int j = 3; 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(int i = 1; i <= m; i++)
        for(int j = 1; j <= m; j++)
        {
            c[i][j] = a[i][j] * (1 - a[i - 1][j]);
            d[i][j] = b[i][j] * (1 - a[i - 1][j]);
        }
    dp[1][1] = 1; dp[1][2] = 2;
    for(int i = 2; i <= m; i++)
    {
        lf sum = 0;
        for(int k = 2; k <= m; k++)
        { 
            dp[i][1] += dp[i - 1][k] * d[i - 1][k];
            sum += d[i - 1][k];
        }
        dp[i][1] = 1 + dp[i][1] / sum;
        for(int j = 2; j <= 50; j++)
        {
            sum = 0;
            for(int k = 1; k < j; k++)
            {
                dp[i][j] += dp[i - 1][k] * c[i - 1][k];
                sum += c[i - 1][k];
            }
            dp[i][j] = j + dp[i][j] / sum;
        }
    }
    if(n <= m)
    {
        lf res = 0;
        for(int i = 1; i <= n + 1; i++) res += dp[n][i] * c[n][i];
        cout << fixed << setprecision(12) << res << '\n'; return 0;
    }
    static mat A(5, vector<lf>(m + 5)), B(m + 5, vector<lf>(m + 5));
    A[0][0] = B[0][0] = 1;
    for(int i = 1; i <= m; i++) A[0][i] = dp[m][i];
    lf sum = 0;
    for(int j = 2; j <= m; j++) { B[j][1] += d[m][j], sum += d[m][j]; }
    for(int j = 2; j <= m; j++) B[j][1] /= sum;
    B[0][1] = 1;
    for(int i = 2; i <= m; i++)
    {
        sum = 0;
        for(int j = 1; j < i; j++) { B[j][i] += c[m][j], sum += c[m][j]; }
        for(int j = 1; j < i; j++) B[j][i] /= sum;
        B[0][i] = i;
    }
    A = mul(A, B = qpow(B, n - m));
    lf res = 0;
    for(int i = 1; i <= m; i++) res += A[0][i] * c[m][i];
    cout << fixed << setprecision(12) << res << '\n';
    return 0;
}