CF1097 Hello2019 D Makoto and a Blackboard 期望 dp

题目链接

题意：
给你一个数$n$，这个数每次会等概率地变成它的任意一个因数，求$k$次操作后期望得到的数字。$n<=1e15,k<=10000$

题解：
先吐槽两句。熬夜打Hello2019，结果自闭。全机房就我rating最低，结果还掉分了。C就是个水题，结果我还脑子一抽WA了3次。这个题被EthanZyh大佬当场用半小时秒掉，然而我想了一个多小时不会。并且更惨的是，想着想着D不知不觉感觉头晕，以为是思考问题导致的，结果第二天发烧到38度。。。

感觉这个题还不错，赛后看了看题解，参考了一下zyh大佬的代码。

这个题首先我们对$n$进行质因数分解，复杂度$O(\sqrt{n})$即可。因为这个题我们没法考虑枚举因数，因数是$\sqrt{n}$量级的，我当时想了很久，发现复杂度都不对。当时有想枚举质因子来做，但是脑子比较晕，没有想清楚。我们来考虑假如$n$是一个质数的若干次方的形式，也就是$n=p^x$的情况，那么我们设$f[i][j]$表示$i$次操作之后$p$的指数是$j$的概率。转移并不难想，我们有$dp[i][j]={\sum }_{k=j}^x\frac{dp[i-1][k]}{k+1}$，除以$k+1$的原因是指数可以随机变为$0-k$的任意一个，有$k+1$种可能。我当时其实就是没有想到对于每一个质因子，你可以独立处理出$k$次操作之后变成每一个指数的概率，然后每一个质因子之间是相互独立的。所以对于一个一般的$n$，我们记录下它有哪些质因子以及每一个质因子有多少个，然后对于每一个质因子去跑一遍dp。接下来就是通过这些质因子之间任意相乘组合出所有$n$的因数，同时你还要知道每一个质因数出现了多少次，这样就可以乘上你之前算的概率了。写法是用一个dfs，有三个状态，第一个是记录当前搜索到第几个质因数，第二个是记录当前以及算过的质因数的若干次方相乘之后暂时得到的数字是多少，第三个是记录$k$次操作后只考虑前面的那些质因数，出现现在这个数字的概率是多少。我们每次枚举当前因数选多少个，进入下一层的时候更新当前数字，并且用dp出来的结果更新当前数字的概率即可。最后期望就是每一个因数出现的概率乘数值之和。

我们设$n$的质因子个数是$x$，那么复杂度$O(\sqrt{n}+k\ast logn\ast x^3)$，虽然看上去后半部分有点不靠谱，但是仔细分析一下，$x$是严格小于$logn$的，并且你如果因子种类多了的话，他对应的次数就会少，所以可能是乘$x^2$的？其实那个应该的$log1e9+7$，再加上$k$只有10000，所以反正是能跑过的。

代码：

#include <bits/stdc++.h>
using namespace std;

long long n,f[10010][400],p[400],dp[400][400],ans;
const long long mod=1e9+7; 
int k,m,num[400];
inline long long ksm(long long x,long long y)
{
	long long res=1;
	while(y)
	{
		if(y&1)
		res=res*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return res;
}
inline void calc(int x)
{
	for(int i=0;i<=k;++i)
	{
		for(int j=0;j<=num[x];++j)
		f[i][j]=0;
	}
	f[0][num[x]]=1;
	for(int i=1;i<=k;++i)
	{
		for(int j=0;j<=num[x];++j)
		{
			long long ji=ksm(j+1,mod-2);
			for(int l=0;l<=j;++l)
			f[i][l]=(f[i][l]+f[i-1][j]*ji%mod)%mod;
		}
	}
}
inline void dfs(int x,long long y,long long z)
{
	if(x>m)
	{
		ans=(ans+y%mod*z%mod)%mod;
		return;
	}
	long long ji=1;
	for(int i=0;i<=num[x];++i)
	{
		dfs(x+1,y*ji,z*dp[x][i]%mod);
		ji*=p[x];
	}
}
int main()
{
	scanf("%I64d%d",&n,&k);
	for(long long i=2;i*i<=n;++i)
	{
		if(n%i==0)
		{
			p[++m]=i;
			while(n%i==0)
			{
				n/=i;
				++num[m];
			}
		}
	}
	if(n>1)
	{
		p[++m]=n;
		num[m]=1;
	}
	for(int i=1;i<=m;++i)
	{
		calc(i);
		for(int j=0;j<=num[i];++j)
		dp[i][j]=f[k][j];
	}
	dfs(1,1,1);
	printf("%I64d\n",ans);
	return 0;
}