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