Wunder Fund Round 2016 (Div. 1 + Div. 2 combined) G

题面：

有一个初始为空的序列，在序列末尾随机添加1或2，有p的概率添加1，1 - p的概率添加2，如果序列末尾有连着的两个相同的数k，那么他们会合并成k + 1，这个合并只要可行，可以一直持续下去

序列有个长度限制n，当序列凑到n位且无法合并，添加操作就结束，问结束时所有权值的和的期望。

solution:

假设数字i出现的概率为p(i)，则数字i + 1出现的期望可以粗略地估计为p(i)^2，那么，值大于50的数出现的概率就非常小了，对答案的贡献可以忽略不计

定义a[i][j]:在长度限制为i的情况下，最左端出现过数字j，最后让序列长度停止在i的概率

那么有a[i][j] = a[i][j - 1] * a[i - 1][j - 1]，即先出现一次j - 1，再出现j - 1的概率

类似地，定义b[i][j]，但加一个约束，序列初始时的第一个数字必须是2，有b[i][j] = b[i][j - 1] * a[i - 1][j - 1]

边界为a[i][1] = p，a[i][2] = b[i][2] = 1 - p

定义f[i][j]:构造一个长度为i的序列，构造结束时序列最左端为j，序列的和的期望

一般地，有f[i][j] = j + ∑(f[i - 1][k] * a[i - 1][k] * (1 - a[i - 2][k])) / ∑(a[i - 1][k] * (1 - a[i - 2][k])) (k = 1 ~ j - 1)

因为第一位确定为j时，第二位出现过的数字一定不能大于等于j，否则会和这个j合并

对于第三位，因为第二位的值已经确定，所以只需保证第三位不为k即可，乘上概率转移一下

特别地，当j = 1，

f[i][1] = 1 + ∑(f[i - 1][k] * b[i - 1][k] * (1 - a[i - 2][k])) / ∑(b[i - 1][k] * (1 - a[i - 2][k])) (k = 2 ~ 50)

最后统计答案为Ans = ∑f[n][i] * a[n][i] * (1 - a[n][i])

注意到转移超过五十次时，a，b数组的值就不在更新了

所以暴力前50次转移，后面的转移就可以用矩阵乘法了

#include<iostream>#include<cstdio>#include<algorithm>#include<cmath>#include<cstring>#include<vector>#include<queue>#include<set>#include<map>#include<stack>#include<bitset>#include<ext/pb_ds/priority_queue.hpp>using namespace std;const int N = 51;typedef long double DB;const DB One = 1.000;struct data{DB a[N][N]; data(){memset(a,0,sizeof(a));}data operator * (const data &b){data c;for (int k = 0; k < N; k++)for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)c.a[i][j] += a[i][k] * b.a[k][j];return c;}};int n;DB Ans,p,q,a[N][N],b[N][N],f[N][N];void Calc_50(){for (int i = 2; i < N; i++){for (int j = 1; j < N; j++){if (j == 1) a[i][j] += p;if (j == 2) a[i][j] += q,b[i][j] += q;a[i][j] += a[i][j - 1] * a[i - 1][j - 1];b[i][j] += b[i][j - 1] * a[i - 1][j - 1];}DB s1,s2; s1 = s2 = 0;for (int k = 2; k <= 50; k++){s2 += b[i - 1][k] * (One - a[i - 2][k]);s1 += f[i - 1][k] * b[i - 1][k] * (One - a[i - 2][k]);}f[i][1] = One + s1 / s2;for (int j = 2; j <= 50; j++){s1 = s2 = 0;for (int k = 1; k < j; k++){s2 += a[i - 1][k] * (One - a[i - 2][k]);s1 += f[i - 1][k] * a[i - 1][k] * (One - a[i - 2][k]);}f[i][j] = j + s1 / s2;}}}data ksm(data x){data ret;for (int i = 0; i < N; i++) ret.a[i][i] = 1;for (; n; n >>= 1){if (n & 1) ret = ret * x;x = x * x;}return ret;}void Calc(){data k; DB sum = 0;f[50][0] = k.a[0][0] = 1;for (int i = 2; i <= 50; i++)sum += b[50][i] * (One - a[50][i]);k.a[0][1] = 1;for (int i = 2; i <= 50; i++)k.a[i][1] = b[50][i] * (One - a[50][i]) / sum;for (int i = 2; i <= 50; i++){sum = 0; k.a[0][i] = (DB)(i);for (int j = 1; j < i; j++) sum += a[50][j] * (One - a[50][j]);for (int j = 1; j < i; j++) k.a[j][i] = a[50][j] * (One - a[50][j]) / sum;}k = ksm(k);for (int i = 1; i <= 50; i++)for (int j = 0; j <= 50; j++)Ans += f[50][j] * k.a[j][i] * a[50][i] * (One - a[50][i]);printf("%.12lf\n",(double)(Ans));}int main(){#ifdef DMCfreopen("DMC.txt","r",stdin);#endifcin >> n >> p; p /= 1E9; q = One - p;a[1][1] = p; a[1][2] = b[1][2] = q;f[1][1] = 1; f[1][2] = 2; Calc_50();if (n <= 50){for (int i = 1; i <= 50; i++)Ans += f[n][i] * a[n][i] * (One - a[n - 1][i]);printf("%.12lf\n",(double)(Ans)); return 0;}n -= 50; Calc();return 0;}