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

最新推荐文章于 2021-01-31 22:48:56 发布

CRZbulabula

最新推荐文章于 2021-01-31 22:48:56 发布

阅读量341

点赞数

CC 4.0 BY-SA版权

分类专栏： dp 矩阵乘法

本文链接：https://blog.youkuaiyun.com/CRZbulabula/article/details/69177524

dp 同时被 2 个专栏收录

152 篇文章

订阅专栏

矩阵乘法

18 篇文章

订阅专栏

本文探讨了一个关于序列构造的问题，通过随机添加1或2并考虑特定的合并规则，求解序列长度达到指定值时所有数值之和的期望。采用概率论与动态规划相结合的方法，实现了高效准确的解答。

题面：

有一个初始为空的序列，在序列末尾随机添加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 DMC
		freopen("DMC.txt","r",stdin);
	#endif
	
	cin >> 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;
}