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