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 _{k=2}^n{f}_{i-1,k}{B}_{i-1,k}}{\sum _{k=2}^m{B}_{i-1,k}}(j=1)$
$f_{i,j}=j+\frac{\sum _{k=1}^{j-1}{f}_{i-1,k}{A}_{i-1,k}}{\sum _{k=1}^{j-1}{A}_{i-1,k}}(j\ne 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⩾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 _{k=2}^{\min (n,m)}{f}_{i-1,k}{B}_{\min (i-1,m),k}}{\sum _{k=2}^{\min (n,m)}{B}_{\min (i-1,m),k}}(j=1)$
$f_{i,j}=j+\frac{\sum _{k=1}^{\min (j-1,m)}{f}_{i-1,k}{A}_{\min (i-1,m),k}}{\sum _{k=1}^{\min (j-1,m)}{A}_{\min (i-1,m),k}}(j\ne 1)$

然后我们如果把 $⩽m$ 的单独暴力处理，然后看 $>m$ 的部分：
$f_{i,j}=j+\frac{\sum _{k=2}^m{f}_{i-1,k}{B}_{m,k}}{\sum _{k=2}^m{B}_{m,k}}(j=1)$
$f_{i,j}=j+\frac{\sum _{k=1}^m{f}_{i-1,k}{A}_{m,k}}{\sum _{k=1}^m{A}_{m,k}}(j\ne 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;
}