【luogu CF618G】Combining Slimes（矩阵乘法）（DP）
Combining Slimes
题目链接：luogu CF618G
题目大意
有一个长度为 n 的栈，如果栈顶两个值都是 x 就会合并成 x+1，一开始没有东西。
你有 p 的概率放进去一个 1，1-p 的概率放入 2，问你当栈被放满的时候，你的期望分数。
分数是栈里所有值的和。
思路
考虑统计每个位置的数的贡献。
由于某个位置的数会因为合并改变，考虑先搞点简单的，某个位置出现过：
\(a_{i,j}\) 为用了 \(i\) 个位置，最左边是 \(j\) 的概率：
初始：\(a_{i,1}=p,a_{i,2}=1-p\)
转移：\(a_{i,j}=a_{i,j-1}a_{i-1,j-1}\)（就是后面的两个位置都是 \(j-1\)，合并成 \(j\)）
发现 \(i>1\) 的时候 \(a_{i,1},a_{i,2}\) 又要初始又要转移，那就特判一下都加上。
那我们要的不是出现过，而是它不会被合并掉，最后留了下来。
\(A_{i,j}\) 为用了 \(i\) 个位置，最左边是 \(j\) 而且不会被合并掉的概率：
初始化：\(A_{1,1}=p,A_{1,2}=1-p\)
\(A_{i,j}=a_{i,j}(1-a_{i-1,j})\)
那我们就可以试着算算期望了：
\(f_{i,j}\) 为最后序列右边 \(i\) 位最左边是 \(j\) 的情况下，这 \(i\) 位的期望和。
那就直接枚举右边 \(i-1\) 位的，考虑一下范围。
不难想想，只有两种可能，要么是你自己是 \(1\)，它不是 \(1\)，要么它就比你小。
\(1\) 那个显然，如果你不是 \(1\)，它比你打一定要从至少 \(2\) 往上一步一步加，那 \(2\) 到它那个数之间的数都要经历过。
那到你那个数的时候就合并了啊！
然后发现 \(2\) 这个好像很特别，似乎要再 DP 一下：
\(b_{i,j}\) 为长度为 \(i\) 的，第一次放进来的是 \(2\)，最左边是 \(j\) 的概率：
初始化：\(b_{i,2}=1-p\)
\(b_{i,j}=b_{i,j-1}a_{i-1,j-1}\)
同样的道理，我们也要求不变的：
\(B_{i,j}\) 为长度为 \(i\)，第一次放进来是 \(2\)，最左边是 \(j\) 而且不会被合并的概率：
初始化：\(B_{1,2}=1-p\)
\(B_{i,j}=b_{i,j}(1-a_{i-1,j})\)
这些都弄好之后，我们正式开始搞 \(f\) 的转移：
\(f_{i,j}=j+\frac{\sum\limits_{k=2}^{n}f_{i-1,k}B_{i-1,k}}{\sum\limits_{k=2}^mB_{i-1,k}}(j=1)\)
\(f_{i,j}=j+\frac{\sum\limits_{k=1}^{j-1}f_{i-1,k}A_{i-1,k}}{\sum\limits_{k=1}^{j-1}A_{i-1,k}}(j\neq 1)\)
那我们就可以转移啦！
但是 \(n\) 很大，那你这个是 \(n^2+nm\) 的啊。
发现一个事情，你会觉得它一直没有 \(1\) 接着 \(2\) 的情况其实会很少。
而且你要注意到你凑出 \(x\) 你至少需要 \(2^{x-2}\) 次。
那比如一个 \(50\)，那不出现 \(1,2\) 的概率就是 \((1-p(1-p))^{2^{50}}\)，只要下面不是 \(1\) 基本上就宣告约等于 \(0\) 了。
那合并的概率也就是 \(0\)。
也就是说，我们设 \(m=50\)，那 \(j\geqslant m\) 的时候，我们可以认为 \(A_{i,j}=0,B_{i,j}=0\)。
然后发现 \(i,j\) 太大的 \(i\) 也不好搞啊，但是你看看那个式子：
\(A_{i,j}=a_{i,j}(1-a_{i-1,j})\)，\(B_{i,j}\) 也差不多。
那 \(j\) 很大的时候，\(a_{i-1,j},a_{i,j}\) 都很小，那 \(A_{i,j}\) 不也很小，那就 \(0\) 咯，直接忽略，所以也只用看到 \(m\)。
（毕竟这题是有精度要求的，不是取模那些）
那再看式子：
\(f_{i,j}=j+\frac{\sum\limits_{k=2}^{\min(n,m)}f_{i-1,k}B_{\min(i-1,m),k}}{\sum\limits_{k=2}^{\min(n,m)}B_{\min(i-1,m),k}}(j=1)\)
\(f_{i,j}=j+\frac{\sum\limits_{k=1}^{\min(j-1,m)}f_{i-1,k}A_{\min(i-1,m),k}}{\sum\limits_{k=1}^{\min(j-1,m)}A_{\min(i-1,m),k}}(j\neq 1)\)
然后我们如果把 \(\leqslant m\) 的单独暴力处理，然后看 \(>m\) 的部分：
\(f_{i,j}=j+\frac{\sum\limits_{k=2}^{m}f_{i-1,k}B_{m,k}}{\sum\limits_{k=2}^{m}B_{m,k}}(j=1)\)
\(f_{i,j}=j+\frac{\sum\limits_{k=1}^{m}f_{i-1,k}A_{m,k}}{\sum\limits_{k=1}^{m}A_{m,k}}(j\neq 1)\)
怎么感觉，有点像那种递推的式子。
看看，会发现确实可以用 \(1\times m\) 的矩阵表示 \(f_i\)。
然后转移矩阵因为 \(i-1\) 都给你变成了 \(j\)，所以是固定的。
那直接上矩阵快速幂就可以啦！
代码
#include<cstdio>
using namespace std;
const int N = 55;
int n, m;
double p, a[N][N], A[N][N], b[N][N], B[N][N], f[N][N];
struct matrix {
int n, m;
double a[N][N];
}A_, B_;
matrix operator *(matrix x, matrix y) {
matrix z; z.n = x.n; z.m = y.m;
for (int i = 0; i < z.n; i++)
for (int j = 0; j < z.m; j++)
z.a[i][j] = 0;
for (int k = 0; k < x.m; k++)
for (int i = 0; i < z.n; i++)
for (int j = 0; j < z.m; j++)
z.a[i][j] += x.a[i][k] * y.a[k][j];
return z;
}
matrix ksm(matrix x, int y) {
matrix z = x; y--;
while (y) {
if (y & 1) z = z * x;
x = x * x; y >>= 1;
}
return z;
}
int main() {
scanf("%d %lf", &n, &p);
p /= 1000000000; m = 50;
a[1][1] = p; a[1][2] = 1 - p;
for (int i = 2; i <= m; i++) {
a[i][1] = p; a[i][2] = 1 - p + a[i][1] * a[i - 1][1];
for (int j = 3; j <= m; j++)
a[i][j] = a[i][j - 1] * a[i - 1][j - 1];
}
A[1][1] = p; A[1][2] = 1 - p;
for (int i = 2; i <= m; i++) {
for (int j = 1; j <= m; j++)
A[i][j] = a[i][j] * (1 - a[i - 1][j]);
}
b[1][2] = 1 - p;
for (int i = 2; i <= m; i++) {
b[i][2] = 1 - p;
for (int j = 3; j <= m; j++)
b[i][j] = b[i][j - 1] * a[i - 1][j - 1];
}
B[1][2] = 1 - p;
for (int i = 2; i <= m; i++) {
for (int j = 1; j <= m; j++)
B[i][j] = b[i][j] * (1 - a[i - 1][j]);
}
f[1][1] = 1; f[1][2] = 2;
for (int i = 2; i <= m; i++) {
double sum = 0;
for (int k = 2; k <= m; k++) {
f[i][1] += f[i - 1][k] * B[i - 1][k]; sum += B[i - 1][k];
}
f[i][1] = f[i][1] / sum + 1;
for (int j = 2; j <= m; j++) {
sum = 0;
for (int k = 1; k < j; k++) {
f[i][j] += f[i - 1][k] * A[i - 1][k];
sum += A[i - 1][k];
}
f[i][j] = f[i][j] / sum + j;
}
}
if (n <= m) {
double ans = 0;
for (int i = 1; i <= m; i++)
ans += A[n][i] * f[n][i];
printf("%.10lf", ans); return 0;
}
A_.n = 1; A_.m = m + 1; A_.a[0][0] = 1;
for (int i = 1; i <= m; i++) A_.a[0][i] = f[m][i];
B_.n = m + 1; B_.m = m + 1; B_.a[0][0] = 1;
for (int i = 1; i <= m; i++) {
B_.a[0][i] = i; if (i == 1) continue; double sum = 0;
for (int j = 1; j < i; j++) sum += A[m][j];
for (int j = 1; j < i; j++) B_.a[j][i] = A[m][j] / sum;
}
double sum = 0;
for (int i = 2; i <= m; i++) sum += B[m][i];
for (int i = 2; i <= m; i++) B_.a[i][1] = B[m][i] / sum;
B_ = ksm(B_, n - m); A_ = A_ * B_;
double ans = 0;
for (int i = 1; i <= m; i++)
ans += A[m][i] * A_.a[0][i];
printf("%.10lf", ans); return 0;
return 0;
}
