【2023NOI模拟赛24】随机除法 做题记录

给定一个正整数 $n$，每次等概率从 $n$ 的所有因子中选择一个，设其为 $d$，将 $n$ 变为 $\frac{n}{d}$。

求期望多少次 $n$ 变为 $1$，对 $10^9+7$ 取模。

$1\le n\le 10^{24}$，输入会给出 $n$ 的所有不同的质因子 $p_{1\sim m}$。

设 $dp_i$ 表示 $n=i$ 时的答案，那么有：

$\begin{aligned} dp_{i}&=\frac{\sum\limits_{d|i,d\le i}dp_d}{k}+1\\ dp_{i}&=\frac{\sum\limits_{d|i,d<i}dp_d}{k}+\frac{dp_i}{k}+1\\ \frac{(k-1)dp_i}{k}&=\frac{\sum\limits_{d|i,d<i}dp_d}{k}+1\\ dp_{i}&=\frac{k+\sum\limits_{d|i,d<i}dp_d}{k-1}\\ \end{aligned}$

设 $n=\prod\limits_{i=1}^m p_i^{e_i}$，那么显然可重集 $\{e_1,e_2,\dots,e_m\}$ 相同的 $n$ 的 $dp_n$ 是相同的。

注意到由于 $\{e_1,e_2,\dots,e_m\}$ 是可重集，所以不妨钦定 $e_i\le e_{i+1}$，不同的 $\{e_1,e_2,\dots,e_m\}$ 最多只有 $\prod\limits_{i=1}^m\log_{\prod\limits_{j=1}^ip_j}(n)$ 个，实际上只有大概 $2\times 10^5$ 个。

那么设 $dp_{\{e_1,e_2,\dots,e_m\}}$ 为可重集为 $\{e_1,e_2,\dots,e_m\}$ 时的答案，则转移是一个高维前缀和，可以设 $f_{\{e_1,e_2,\dots,e_m\},i}=\sum\limits_{\{e'_1,e'_2,\dots,e'_m\}}\prod\limits_{j=1}^i[e'_j\le e_j]\prod\limits_{j=i+1}^m[e'_j=e_j]dp_{\{e'_1,e'_2,\dots,e'_m\}}$，则有转移：

$f_{\{e_1,e_2,\dots,e_m\},i}=f_{\{e_1,e_2,\dots,e_m\},i-1}+f_{\{e_1,e_2,\dots,e_i-1,\dots,e_m\},i}+f_{\{e_1,e_2,\dots,e_m\},0}$

则 $dp_{\{e_1,e_2,\dots,e_m\}}=f_{\{e_1,e_2,\dots,e_m\},0}$。

由于转移需要的是 $f_{\{e_1,e_2,\dots,e_m\},m}-dp_{\{e_1,e_2,\dots,e_m\}}$，所以不妨先令 $f_{\{e_1,e_2,\dots,e_m\},0}=0$，做完高维前缀和再加回去。

时间复杂度 $O(2\times 10^5\times 20)$。

代码如下：

#include <iostream>
#include <algorithm>
#include <cstring>
#include <vector>
#include <map>

using namespace std;

const int SS=105,STS=200005,S=25,p=1000000007;

struct sta
{
	vector<int> v;
	inline void push(int x){v.push_back(x);}
	inline int size(){return v.size();}
	inline void sort(){::sort(v.begin(),v.end());}
	inline int& operator[](int x){return v[x];}
	inline void delbeg(){v.erase(v.begin());}
	inline bool operator<(const sta &b)const
	{
		if(v.size()!=b.v.size()) return v.size()<b.v.size();
		for(int i=0;i<v.size();i++) if(v[i]!=b.v[i]) return v[i]<b.v[i];
		return false;
	}
};

char str[SS];
__int128 n,a[S];
int m;
int tot;
map<sta,int> idx;
int f[STS][S];

inline int qpow(int x,int y)
{
	int res=1;
	for(;y>0;y>>=1,x=1ll*x*x%p) res=y&1?1ll*res*x%p:res;
	return res;
}

int DP(sta &st)
{
	if(idx.find(st)!=idx.end()) return f[idx[st]][0];
	idx[st]=++tot;
	int id=tot;
	if(st.size()==0) return f[id][0]=0;
	for(int i=1;i<=st.size();i++)
	{
		sta nst=st;
		nst[i-1]--,nst.sort();
		if(nst[0]!=0) DP(nst),f[id][i]=(f[id][i-1]+f[idx[nst]][i])%p;
		else nst.delbeg(),DP(nst),f[id][i]=(f[id][i-1]+f[idx[nst]][i-1])%p;
	}
	int ml=1;
	for(int i=0;i<st.size();i++) ml=1ll*ml*(st[i]+1)%p;
	f[id][0]=1ll*qpow(ml-1+p,p-2)*(f[id][st.size()]+ml)%p;
	for(int i=1;i<=st.size();i++) f[id][i]=(f[id][i]+f[id][0])%p;
	return f[id][0];
}

int main()
{
	freopen("div.in","r",stdin);
	freopen("div.out","w",stdout);
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
	while(cin>>str>>m)
	{
		int len=strlen(str);
		n=0;
		for(int i=0;i<len;i++) n=n*10+str[i]-'0';
		sta st;
		for(int i=1;i<=m;i++)
		{
			cin>>str;
			len=strlen(str);
			a[i]=0;
			for(int j=0;j<len;j++) a[i]=a[i]*10+str[j]-'0';
			int cnt=0;
			while(n%a[i]==0) n/=a[i],cnt++;
			st.push(cnt);
		}
		st.sort();
		cout<<DP(st)<<'\n';
	}
	return 0;
}