@codeforces - 618G@ Combining Slimes

---

@description@

一行上摆有 n 个方格。每一次你可以在最右边的方格滴入一滴史莱姆。有 p 的概率该史莱姆大小为 1，有 (1 - p) 的概率该史莱姆大小为 2。
史莱姆会不断往左滚动，直到遇到另一个史莱姆或边界。假如遇到的是大小相同的史莱姆则合并，大小加一，继续往左滚动；否则直接停下。
等到无法操作时，问最后方格内的史莱姆大小总和的期望值。

input：
两个整数 n，p。（1 <= n <= 10^9，1 <= p < 10^9）
其中对应的概率为 p/10^9 与 1 - p/10^9。

output：
输出期望值。

sample input：
2 500000000
sample output：
3.562500000000000
sample explain：
滴入大小为 1 与 2 的史莱姆的概率都为 1/2。
1/4 的概率最终状态为 1 2，3/8 最终为 2 1，3/16 最终为 3 2 与 3/16 最终为 3 1。
期望值为 \(1/4*(1 + 2) + 3/8*(2 + 1) + 3/16*(3 + 2) + 3/16*(3 + 1) = 3.5625\)

@solution@

这道题有点儿意思。

@part - 0@

嘛……首先发现当一个大小为 1 的史莱姆在一个大小为 2 的史莱姆左边，则它们以及它们左边的史莱姆永远不会再变动的。
然后呢？然后我就不会做了。

@part - 1@

根据题解所写，我们要发现这样一个关键的性质：
假如我们需要弄出来大小为 s 的史莱姆，我们需要先弄出来两个大小为 (s-1) 的史莱姆。
令弄出大小为 s 的概率为 \(pr(s)\)，则有递推公式 \(pr(s) = pr^2(s-1)\)，边界 \(pr(2) = (1-p) + p^2\)。
这是一个理想的上界，实际概率可能还要小些。
但是就这种情况下，取 \(pr(2)_{max}=(1-\frac{1}{10^9})\)（不是准确值），则当 s >= 50 时， \(pr(s)\) 已经趋近于无穷小了（可以自行算），不会再对答案产生影响。

所以史莱姆的大小可以控制到 50 以下。

……
MMP 鬼才会去想这个性质。

@part - 2@

其实下面这个性质并不难想到，但是如果没有上面那个性质也运用不了。令 \(pr(s, n)\)表示 n 个方格弄出大小为 s 的史莱姆的概率。
可以发现，当 n 充分大时， \(pr(s, n) = pr(s, n+1)\)。

感性理解一下：你方格数量多了，反而有些方格可能用不到。将这些用不到的方格去掉，概率依然还是不变的。

这个充分大是多大呢？
\(pr(s, n)\) 之间有递推关系 \(pr(s, n) = pr(s-1, n)*pr(s-1, n-1)\)。
边界条件 \(\begin{cases}pr(1, 1) = p, pr(2, 1) = 1-p\\pr(1, n) = p, pr(2,n) = (1-p)+p^2&n\not = 1\end{cases}\)
我们可以发现边界条件其实跟 n 没有关系。然后递推关系中，每一次 s 都要减少 1，所以 s 过后就会碰到边界条件。
所以，当 n > s 时，总有 \(pr(s, n) = pr(s, n+1)\)。

因此只要我们已知 \(pr(s, n)(s \le 50, n \le 50)\) 的值，就可以表示所有的 \(pr(s, n)\) 的值。

然而还有个问题，倒回到我们一开始说的那个性质：有可能我们的最左边固定的是一个 1 2 的形式。
为了解决这个问题，我们令 \(pr'(s, n)\) 表示 n 个方格，最左边已经有一个大小为 2 的史莱姆的前提下，弄出大小为 s 的史莱姆的概率。
类似地可以得到递推关系：\(pr'(s, n)=pr'(s-1,n)*pr(s-1,n-1)\)。
边界条件不再复述，可以自行推导，也可以等会儿看代码。

@part - 3@

发现的这些性质……所以到底有啥子用？
既然上文中有这么多递推关系，我们就来想想用我们的 dp 来搞这道题。

令 \(dp(s, n)\) 表示最终局面从右往左数前 n 个格子，第 n 个格子的史莱姆大小为 s 的前提下，的期望大小和。根据上文， s <= 50。

令 \(f(s, n)=pr(s, n)*(1-pr(s, n-1))\)，含义为第 n 个格子史莱姆大小为 s，第 n-1 个格子史莱姆大小 < s 的概率。这样第 n 个格子史莱姆绝对不会与第 n-1 个合并，只要第 n+1 个格子的史莱姆不与第 n 个格子合并，整个局面就是“稳定的”。

转移时，枚举第 i+1 个格子大小为 k，则：
\[dp(s,n)= s+\sum_{k=1}^{s-1}(dp(k, n-1)*\dfrac{f(k,n-1)}{\sum_{p=1}^{s-1}f(p,n-1)})\]
比较难理解的可能就是右边乘上来的概率为什么长成那样。

其实这（应该）是一个条件概率。
条件概率的公式为 \(P(B|A)=\dfrac{P(AB)}{P(A)}\)，即在事件 A 发生的前提下，事件 B 发生的概率等于事件 A B 同时发生的概率除以 A 发生的概率。
对应到该题。事件 A 为 “最终局面中第 i 个格子的史莱姆大小为 j ”，事件 B 为 “第 i+1 个格子的史莱姆大小为 k”。
嗯。我相信它是这样的。

为什么不直接用 \(dp(s, n)\) 乘上一个 \(f(s, n)\) 呢？因为当 j = 1 时，上面那个转移式，就不能那么转移了。因为大小为 1 后面可能接个 2。
我们令 \(g(s, n)=pr'(s, n)*(1-pr(s, n-1))\)。则：
\[dp(1, n)=1+\sum_{k=1}^{50}(dp(k, n-1)*\dfrac{g(k,n-1)}{\sum_{p=1}^{50}g(p,n-1)})\]

然后，由于我们最后发现的那个性质，所以当 n 充分大时右边乘上来的概率是不会变化的，是个常数。
所以我们就可以跑矩阵加速了。

最后得到 \(ans = \sum_{i=1}^{50}f(i, n)*dp(i, n)\)。

【我才不会说题解基本全靠口胡和感性认知 QAQ】
【因为我自己也不是太懂这道题 QAQ】

@accepted code@

通过浮点数的精度，使得某项变量因为远远超过题目所给的精度而忽略它。
这种类型我还只在 NOI2016 的旷野大计算那道题见过，果然还是太弱 QAQ。

//代码极丑无比，请勿模仿。
#include<cstdio>
const int MAXN = 50;
double pr1[MAXN + 5][MAXN + 5], pr2[MAXN + 5][MAXN + 5];
double f[MAXN + 5][MAXN + 5], g[MAXN + 5][MAXN + 5];
double dp[MAXN + 5][MAXN + 5];
struct matrix{
    double m[MAXN + 5][MAXN + 5];
    int r, c;
}M, R;
void init(double p) {
    pr1[1][1] = p, pr1[2][1] = 1-p;
    f[1][1] = p, f[2][1] = 1-p, g[2][1] = 1;
    for(int i=2;i<=MAXN;i++) {
        pr1[1][i] = p, pr1[2][i] = (1-p) + p*p, pr2[2][i] = 1;
        for(int j=3;j<=MAXN;j++)
            pr1[j][i] = pr1[j-1][i]*pr1[j-1][i-1], pr2[j][i] = pr2[j-1][i]*pr1[j-1][i-1];
        for(int j=1;j<=MAXN;j++)
            f[j][i] = pr1[j][i]*(1-pr1[j][i-1]), g[j][i] = pr2[j][i]*(1-pr1[j][i-1]);
    }
    dp[1][1] = 1, dp[2][1] = 2;
    for(int i=2;i<=MAXN;i++) {
        double del = 0; dp[1][i] = 1;
        for(int l=1;l<=50;l++)
            del += g[l][i-1];
        for(int k=1;k<=MAXN;k++)
            dp[1][i] += dp[k][i-1]*g[k][i-1]/del;
        del = 0;
        for(int j=2;j<=MAXN;j++) {
            dp[j][i] = j; del += f[j-1][i-1];
            for(int k=1;k<j;k++)
                dp[j][i] += dp[k][i-1]*f[k][i-1]/del;
        }
    }
    M.r = M.c = MAXN+1;
    double del = 0;
    for(int i=1;i<=MAXN;i++)
        del += g[i][MAXN];
    for(int i=1;i<=MAXN;i++)
        M.m[1][i] = g[i][MAXN]/del;
    del = 0;
    for(int i=2;i<=MAXN;i++) {
        del += f[i-1][MAXN];
        for(int j=1;j<i;j++)
            M.m[i][j] = f[j][MAXN]/del;
        for(int j=i;j<=MAXN;j++)
            M.m[i][j] = 0;
    }
    for(int i=1;i<=MAXN;i++)
        M.m[0][i] = 0, M.m[i][0] = i;
    M.m[0][0] = 1;
    R.r = MAXN+1, R.c = 1;
    for(int i=1;i<=MAXN;i++)
        R.m[i][0] = dp[i][MAXN];
    R.m[0][0] = 1;
}
matrix operator * (matrix A, matrix B) {
    matrix C; C.r = A.r, C.c = B.c;
    for(int i=0;i<C.r;i++)
        for(int j=0;j<C.c;j++)
            C.m[i][j] = 0;
    for(int i=0;i<A.r;i++)
        for(int j=0;j<B.c;j++)
            for(int k=0;k<A.c;k++)
                C.m[i][j] += A.m[i][k] * B.m[k][j];
    return C;
}
matrix quick_pow(matrix b, int p) {
    matrix ret; ret.r = ret.c = b.r;
    for(int i=0;i<ret.r;i++)
        for(int j=0;j<ret.c;j++)
            ret.m[i][j] = (i == j);
    while( p ) {
        if( p & 1 ) ret = ret * b;
        b = b * b;
        p >>= 1;
    }
    return ret;
}
int main() {
    int n, p;
    scanf("%d%d", &n, &p);
    init(p/1E9);
    if( n <= MAXN ) {
        double ans = 0;
        for(int i=1;i<=MAXN;i++)
            ans += f[i][n]*dp[i][n];
        printf("%lf\n", ans);
    }
    else {
        R = quick_pow(M, n-MAXN)*R;
        double ans = 0;
        for(int i=1;i<=MAXN;i++)
            ans += f[i][MAXN]*R.m[i][0];
        printf("%lf\n", ans);
    }
}

@details@

不知道为什么想这道题的时候脑子里有一堆史莱姆在滚来滚去。好可爱来着 p(#￣▽￣#)o。

我看了好久才发现那个转移式是个条件概率 QAQ。可能有其他的理解方法，但是我太弱了真的想不到，只能用条件概率去解释 QAQ。
各位过路的 dalao 如果能提供更简洁的理解思路麻烦留言在下面好吗 QAQ。
救救蒟蒻吧 QAQ。

转载于:https://www.cnblogs.com/Tiw-Air-OAO/p/10181690.html