【XSY3657】因数分解（容斥，DP）

考虑没有 \(b_i\) 的限制怎么做。

先把 \(n!\) 质因数分解：\(n!=\prod\limits_{i=1}^kp_i^{q_i}\)。

设 \(b_i=\prod\limits_{j=1}^kp_j^{x_{i,j}}\)，那么 \(b_i^{a_i}=\prod\limits_{j=1}^kp_j^{a_ix_{i,j}}\)，\(\prod\limits_{i=1}^m b_i^{a_i}=\prod\limits_{j=1}^kp_j^{\sum\limits_{i=1}^ma_ix_{i,j}}\)。

每一个质因数 \(p_j\) 是独立的，于是我们对于每一个 \(q_j\) 求出把 \(q_j\) 分解成 \(\sum\limits_{i=1}^ma_ix_{i}\) 的方案数再乘起来即可。

求这个方案数是一个背包 DP 的事情，具体实现上我们 DP 求出 \(f_v\) 表示将 \(v\) 分解成 \(\sum\limits_{i=1}^ma_ix_{i}\) 的方案数。

那么现在加上 \(b_i\) 的限制后我们就考虑容斥。

对于某一种 \(\{b_i\}\) 的方案，将所有的 \(i\) 划分成若干组，每一组里的 \(b_i\) 是相同的，但不保证所有相同的 \(b_i\) 划分进同一组，然后设一共有 \(B\) 组，分别是 \(s_1,s_2,\cdots,s_B\)。 那么满足 \(\forall i,j\in s_x,b_i=b_j\)。

考虑设容斥系数 \(f(s_i)\) 使得：

\[ans=\sum_{枚举划分T=(s_1,s_2,\cdots,s_B)}\left(\prod_{i=1}^Bf(|s_i|)\right)dp(T) \]

其中 \(dp(T)\) 表示满足划分 \(T\) 的 \(b\) 的方案数。

换言之，设 \(suma_i=\sum\limits_{j\in s_i}a_j\)，我们设 \(s_i\) 这一组内是都有相同的 \(b'_i\)，此时组已经是分好了的，如果我们能确定 \(b'_i\)，我们就能确定原来所有的 \(b_i\)，也就是我们要求 \({b_i'}\) 使得 \(\prod\limits_{i=1}^B{b_i'}^{suma_i}=n!\)，那么 \(dp(T)\) 就是 \({b_i'}\) 的方案数。

注意这里的 \(b_i'\) 并没有要求两两不同，所以我们仍然可以用背包 DP 来求解。

这里插一小段讲一下怎么求解：由于对于不同的 \(T\) 它的 \(suma_i\) 是有可能不同的，对于每一种 \(T\) 都跑一遍背包 DP 显然不现实，那么我们就不妨直接 dfs 枚举 \(suma\) 数组是什么（设 \(M=\sum a_i\)，那么只需满足 \(\sum_{i} suma_i=M\) 即可），然后一边枚举一遍跑背包 DP，并对于所有枚举的 \(suma\) 数组记录答案。

\(dp(T)\) 求解完成，考虑怎么求容斥系数。

那么考虑一种极大的划分 \(T\)（“极大的” 指这种划分方案中所有相同的 \(b_i\) 都在同一组，即不存在两组的 \(b_i\) 相同），对于某种符合划分 \(T\) 的 \(b\) 的方案 \(S\)，考虑 \(S\) 在上面答案计算式等号左右的贡献：

• 显然这种方案对 \(ans\) 有贡献当且仅当没有一组有多于 \(1\) 个元素（即不存在 \(i\neq j\) 且 \(b_i=b_j\)），故对等号左边的贡献为 \([s_1=s_2=\cdots=s_B=1]\)。

• 等号右边，如果在枚举另一种划分 \(T'\) 时 \(S\) 也被统计进去了（即 \(S\) 也符合划分 \(T'\)），当且仅当能将 \(T\) 的每一组再划分成若干组（可以划分为 \(1\) 组，即不变）而形成划分 \(T'\)。

  于是 \(S\) 在等号右边的贡献为：（定义 \(A=B+C\) 当且仅当 \(A=B\cup C\) 且 \(B\cap C=\emptyset\)，即能看作将集合 \(A\) 划分成两个集合 \(B,C\)）

  \[\left(\sum_{s_1=s_{1,1}+\cdots+s_{1,k_1}}\dfrac{1}{k_1!}\prod_{i=1}^{k_1}f(|s_{1,i}|)\right)\left(\sum_{s_2=s_{2,1}+\cdots+s_{2,k_2}}\dfrac{1}{k_2!}\prod_{i=1}^{k_2}f(|s_{2,i}|)\right)\cdots\left(\sum_{s_B=s_{B,1}+\cdots+s_{B,k_B}}\dfrac{1}{k_B!}\prod_{i=1}^{k_B}f(|s_{B,i}|)\right)\\ \]

  （除以 \(k_1!\) 是因为对于某一种 \(s_1=s_{1,1}+\cdots+s_{1,k_1}\) 它会被计算 \(k_1\) 次）

  发现与具体划分的集合长啥样无关，只和划分的集合的大小有关，于是转为枚举集合的大小：

  \[\prod_{i=1}^B\left(\sum_{|s_i|=x_{i,1}+\cdots+x_{i,k_i}}\dfrac{1}{k_i!}\dbinom{|s_i|}{x_{i,1},\cdots,x_{i,k_i}}\prod_{j=1}^{k_i}f(x_{i,j})\right) \]

  从上一步枚举集合跳到这一步枚举集合大小看似非常直接，但系数仍为 \(\dfrac{1}{k_i!}\) 的原因实际上需要理性证明：

  考虑 \(|s_i|\) 的某种分解：\(|s_i|=w_1y_1+w_2y_2+\cdots+w_cy_c\)，其中 \(y_i\) 互不相同。显然此时 \(k_i=\sum\limits_{i=1}^c w_i\)。

  那么这种情况会被 \(|s_i|=x_{i,1}+\cdots+x_{i,k}\) 枚举 \(\dfrac{k_i!}{w_1!w_2!\cdots w_c!}\) 次，那么我们就需要除掉 \(\dfrac{k_i!}{w_1!w_2!\cdots w_c!}\)。

  同时，这种情况在每次被枚举到时同一种分配方案都会因为 \(\dbinom{|s_i|}{x_{i,1},\cdots,x_{i,k_i}}\) 而被算重 \(w_1!w_2!\cdots w_c!\) 次，那么我们就需要除掉 \(w_1!w_2!\cdots w_c!\)。

  于是需要除掉 \(\dfrac{k_i!}{w_1!w_2!\cdots w_c!}\times w_1!w_2!\cdots w_c!=k_i!\) 。

  化简一下：

  \[\begin{aligned} &\prod_{i=1}^B\left(\sum_{|s_i|=x_{i,1}+\cdots+x_{i,k_i}}\dfrac{1}{k_i!}\dbinom{|s_i|}{x_{i,1},\cdots,x_{i,k_i}}\prod_{j=1}^{k_i}f(x_{i,j})\right)\\ =&\prod_{i=1}^B\left(\sum_{|s_i|=x_{i,1}+\cdots+x_{i,k_i}}\dfrac{|s_i|!}{k_i!}\prod\limits_{j=1}^{k_i}\dfrac{f(x_{i,j})}{x_{i,j}!}\right)\\ \end{aligned} \]

现在令等号左右贡献相等：

\[[|s_1|=|s_2|=\cdots |s_B|=1]=\prod_{i=1}^B\left(\sum_{|s_i|=x_{i,1}+\cdots+x_{i,k_i}}\dfrac{|s_i|!}{k_i!}\prod\limits_{j=1}^{k_i}\dfrac{f(x_{i,j})}{x_{i,j}!}\right) \]

显然如果下式满足，上式也满足：

\[[|s|=1]=\sum_{|s|=x_1+x_2+\cdots x_{k}}\dfrac{|s|!}{k!}\prod_{j=1}^k\dfrac{f(x_j)}{x_j!} \]

注意到这里有个 EGF 的形式，即如果我们设 \(F(x)=\sum\limits_{i\geq 1}\dfrac{f(i)}{i!}x^i\)，那么就有：

\[\begin{aligned} \big[|s|=1\big]=&\sum_{|s|=x_1+x_2+\cdots x_{k}}\dfrac{|s|!}{k!}\prod_{j=1}^k\dfrac{f(x_j)}{x_j!}\\ =&|s|!\sum_{k\geq 1}\dfrac{1}{k!}[x^{|s|}]F^k(x)\\ =&|s|![x^{|s|}]\left(\sum_{k\geq 0}\dfrac{1}{k!}F^k(x)-1\right)\\ =&|s|![x^{|s|}]\left(e^{F(x)}-1\right)\\ \end{aligned} \]

于是：

\[\begin{aligned} \big[x^{|s|}\big]\left(e^{F(x)}-1\right)&=\dfrac{[|s|=1]}{|s|!}\\ e^{F(x)}&=x+1\\ F(x)&=\ln(x+1) \end{aligned} \]

考虑通过求导拆掉 \(\ln\)：

\[\begin{aligned} F'(x)&=\dfrac{1}{x+1}=\sum_{i\geq 0}(-1)^{i}x^i\\ F(x)&=\sum_{i\geq 1}\dfrac{(-1)^{i-1}}{i}x^i \end{aligned} \]

根据 \(F(x)\) 的定义就有：

\[\dfrac{f(i)}{i!}=\dfrac{(-1)^{i-1}}{i} \]

得到 \(f(i)=(-1)^{i-1}(i-1)!\)。

代入原来的式子求 \(ans\) 即可。

至于 DP 的具体过程，注意到我们每次枚举需要的只是所有 \(s_i\) 的大小和所有的 \(suma_i\)，于是可以以这两个东西当做维度来 DP。

具体过程详见代码：

#include<bits/stdc++.h>

#define A 35
#define N 10010

using namespace std;

namespace modular
{
	const int mod=1000000007;
	inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
	inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
	inline int mul(int x,int y){return 1ll*x*y%mod;}
	inline void Add(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
	inline void Mul(int &x,int y){x=1ll*x*y%mod;}
}using namespace modular;

inline int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

struct data
{
	int size,sum;
	data(){};
	data(int a,int b){size=a,sum=b;}
};

bool operator < (data a,data b)
{
	if(a.sum==b.sum) return a.size<b.size;
	return a.sum<b.sum;
}

int n,m,M,a[A];
int cnt,prime[N];
int maxq,q[N],numq[N];
bool notprime[N];
int fac[A],ifac[A];

map<vector<int>,int>f1;//f1[vec]：当suma数组为vec时，有多少种b的方案（b可以相同）
map<vector<data>,int>f2[2];//f2[pair(size,sum)]：当|si|数组为size、suma数组为sum时，划分T的方案数

int calc(int p)
{
	int ans=0;
	int x=n;
	while(x)
	{
		ans+=x/p;
		x/=p;
	}
	return ans;
}

void init()
{
	fac[0]=1;
	for(int i=1;i<=m;i++) fac[i]=mul(fac[i-1],i);
	ifac[m]=poww(fac[m],mod-2);
	for(int i=m;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
	for(int i=2;i<=n;i++)
	{
		if(!notprime[i])
		{
			prime[++cnt]=i;
			q[i]=calc(i);
			numq[q[i]]++;
			maxq=max(maxq,q[i]);
		}
		for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
		{
			notprime[i*prime[j]]=1;
			if(!(i%prime[j])) break;
		}
	}
}

int aa[A];
int dp[A][N];

void dfs(int k,int lst,int sum)
{
	if(sum==M)
	{
		int ans=1;
		for(int i=1;i<=maxq;i++)
			Mul(ans,poww(dp[k-1][i],numq[i]));
		vector<int>now;
		for(int i=1;i<k;i++)
			now.push_back(aa[i]);
		f1[now]=ans;
		return;
	}
	for(int i=lst;sum+i<=M;i++)
	{
		aa[k]=i;
		for(int j=0;j<=maxq;j++) dp[k][j]=dp[k-1][j];
		for(int j=i;j<=maxq;j++) Add(dp[k][j],dp[k][j-i]);
		dfs(k+1,i,sum+i);
	}
}

int main()
{
	n=read(),m=read();
	init();
	for(int i=1;i<=m;i++)
	{
		a[i]=read();
		M+=a[i];
	}
	dp[0][0]=1;
	dfs(1,1,0);
	bool lst=0,now=1;
	vector<data>empty;
	f2[now][empty]=1;
	for(int i=1;i<=m;i++)
	{
		swap(now,lst);
		f2[now].clear();
		for(auto it=f2[lst].begin();it!=f2[lst].end();it++)
		{
			vector<data>T=it->first;
			for(int j=0;j<(int)T.size();j++)
			{
				vector<data>nT=T;
				nT[j].sum+=a[i],nT[j].size++;
				sort(nT.begin(),nT.end());
				Add(f2[now][nT],it->second);
			}
			vector<data>nT=T;
			nT.push_back(data(1,a[i]));
			sort(nT.begin(),nT.end());
			Add(f2[now][nT],it->second);
		}
	}
	int ans=0;
	for(auto it=f2[now].begin();it!=f2[now].end();it++)
	{
		vector<data>T=it->first;
		int coef=1;
		for(data s:T) Mul(coef,mul((s.size&1)?1:mod-1,fac[s.size-1]));
		vector<int>tmp;
		for(data s:T) tmp.push_back(s.sum);
		Add(ans,mul(coef,mul(f1[tmp],it->second)));
	}
	sort(a+1,a+m+1);
	int tmp=1;
	for(int i=2;i<=m;i++)
	{
		if(a[i]==a[i-1]) tmp++;
		else
		{
			ans=mul(ans,ifac[tmp]);
			tmp=1;
		}
	}
	ans=mul(ans,ifac[tmp]);
	printf("%d\n",ans);
	return 0;
}
/*
10 6
1 2 2 3 3 3
*/
/*
6 4
1 1 2 1
*/