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

考虑没有 $b_i$ 的限制怎么做。

先把 $n!$ 质因数分解：$n!=\prod _{i=1}^k{p}_i^{q_i}$。

设 $b_i=\prod _{j=1}^k{p}_j^{x_{i,j}}$，那么 $b_i^{a_i}=\prod _{j=1}^k{p}_j^{a_i{x}_{i,j}}$，$\prod _{i=1}^m{b}_i^{a_i}=\prod _{j=1}^k{p}_j^{\sum _{i=1}^m{a}_i{x}_{i,j}}$。

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

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

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

对于某一种 { b i } \{b_i\} {bi​} 的方案，将所有的 $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}^{B}f(|s_i|)\right)dp(T)$$
其中 $dp(T)$ 表示满足划分 $T$ 的 $b$ 的方案数。

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

注意这里的 $b_i^{\prime }$ 并没有要求两两不同，所以我们仍然可以用背包 DP 来求解。

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

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

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

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

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

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

  $$\left(\sum _{s_1=s_{1,1}+\cdots +s_{1,k_1}}\frac{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}}\frac{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}}\frac{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}}\frac{1}{k_i!}\left(\genfrac{}{}{0ex}{}{|s_i|}{x_{i,1},\cdots ,x_{i,k_i}}\right)\prod _{j=1}^{k_i}f(x_{i,j})\right)$$
  从上一步枚举集合跳到这一步枚举集合大小看似非常直接，但系数仍为 $\frac{1}{k_i!}$ 的原因实际上需要理性证明：

  考虑 $|s_i|$ 的某种分解：$|s_i|=w_1{y}_1+w_2{y}_2+\cdots +w_c{y}_c$，其中 $y_i$ 互不相同。显然此时 $k_i=\sum _{i=1}^c{w}_i$。

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

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

  于是需要除掉 $\frac{k_i!}{w_1!w_2!\cdots w_c!}\times w_1!w_2!\cdots w_c!=k_i!$ 。

  化简一下：
  $$\begin{array}{rl}& \prod _{i=1}^B\left(\sum _{|s_i|=x_{i,1}+\cdots +x_{i,k_i}}\frac{1}{k_i!}\left(\genfrac{}{}{0ex}{}{|s_i|}{x_{i,1},\cdots ,x_{i,k_i}}\right)\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}}\frac{|s_i|!}{k_i!}\prod _{j=1}^{k_i}\frac{f(x_{i,j})}{x_{i,j}!}\right)\end{array}$$

现在令等号左右贡献相等：
$$[|s_1|=|s_2|=\cdots |s_B|=1]=\prod _{i=1}^B\left(\sum _{|s_i|=x_{i,1}+\cdots +x_{i,k_i}}\frac{|s_i|!}{k_i!}\prod _{j=1}^{k_i}\frac{f(x_{i,j})}{x_{i,j}!}\right)$$
显然如果下式满足，上式也满足：
$$[|s|=1]=\sum _{|s|=x_1+x_2+\cdots x_k}\frac{|s|!}{k!}\prod _{j=1}^k\frac{f(x_j)}{x_j!}$$
注意到这里有个 EGF 的形式，即如果我们设 $F(x)=\sum _{i\ge 1}\frac{f(i)}{i!}{x}^i$，那么就有：
$$\begin{array}{rl}[|s|=1]=& \sum _{|s|=x_1+x_2+\cdots x_k}\frac{|s|!}{k!}\prod _{j=1}^k\frac{f(x_j)}{x_j!}\\ =& |s|!\sum _{k\ge 1}\frac{1}{k!}[x^{|s|}]F^k(x)\\ =& |s|![x^{|s|}]\left(\sum _{k\ge 0}\frac{1}{k!}{F}^k(x)-1\right)\\ =& |s|![x^{|s|}]\left(e^{F(x)}-1\right)\end{array}$$
于是：
$$\begin{array}{rl}\left[x^{|s|}\right]\left(e^{F(x)}-1\right)& =\frac{[|s|=1]}{|s|!}\\ e^{F(x)}& =x+1\\ F(x)& =\ln (x+1)\end{array}$$
考虑通过求导拆掉 $\ln$：
$$\begin{array}{rl}{F}^{\prime }(x)& =\frac{1}{x+1}=\sum _{i\ge 0}(-1{)}^i{x}^i\\ F(x)& =\sum _{i\ge 1}\frac{(-1{)}^{i-1}}{i}{x}^i\end{array}$$
根据 $F(x)$ 的定义就有：
$$\frac{f(i)}{i!}=\frac{(-1{)}^{i-1}}{i}$$
得到 $f(i)=(-1{)}^{i-1}(i-1)!$。

代入原来的式子求 $ans$ 即可。

至于 DP 的具体过程，注意到我们每次枚举需要的只是所有 $s_i$ 的大小和所有的 $sum{a}_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
*/