【XSY3657】 因数分解
DP+容斥神仙题！！！
(b_i)要互不相同，这很烦，于是我们就先考虑让(b_i)可以相同，再容斥。
另外，要让每个元素都不相同，最后再消去顺序。
首先我们考虑到(sum m)很小，于是我们可以先对(n!)进行质因数分解，然后用背包DP算出每一种(sum b_i=m)的划分的分配质因数的方案数。具体来说，我们对每一种(b_i)的划分，考虑把每一个(b_i)当成元素且有无限个，然后进行背包DP，每一个次数为(q_i)的(p_i)的贡献就会是由(b_i)组成(q_i)的方案数。
之后我们考虑对于题目给定的(a_i)，把其中钦定(b_i)相同的划分到同一组，并统计一下方案数，这个也可以用DP来解决。
那么，对于每一种方案，它的贡献就是(分配质因数的方案数 imes 划分a_i的方案数 imes 容斥系数)。
于是来考虑容斥系数。
不妨大胆假设，容斥系数只和划分到同一组的(a_i)的个数有关系，为(f(a_i))。
我们考虑对于每一种划分，我们只钦定了同一组里面的(b_i)相同，并未钦定不同组的(b_i)不同，于是考虑对于每一个(b_i)互不相同的划分，每一个(b_i)相同的组的大小为(s_i)，那么这个划分的贡献有：
[[s_1=s_2=...=s_k=1]=sum_{s_i划分为s'_1,s'_2,...,s'_{k_i}} prod_{i=1}^k frac{s_i!prod_{j=1}^{k_i}f(s'_j)}{k_i!prod_{j=1}^{k_i}s'_j!}
]
其中，“(s_i)被划分为...”就意味着全部统计了(s_i)这种(b_i)互不相同划分的方案数。
右边的那个阶乘式的贡献还是很显然的。
考虑我们把每一个(s)提出来：
[[s=1]=sum_{s划分为s_1,s_2,...,s_k} frac{s!prod_{i=1}^kf(s_i)}{k!prod_{i=1}^{k}s_i!}
]
看到右边的右半部分，不难想到这是一个(EGF)的形式，于是设(F(x))为(f)的(EGF)，则有：
[ecause F(x)=sum_{i=0}^{infty} frac{f(i)}{i!}x^i\
herefore e^{F(x)}=sum_{i=0}^{infty} frac{F^i(x)}{i!}=sum_{i=0}^{infty}frac{(sum_{j=0}^{infty} frac{f(j)}{j!})^i}{i!}\
herefore [x^m]e^{F(x)}=sum_{x_1+x_2+...+x_k=m} frac{1}{k!}prod_{i=1}^kfrac{f(x_i)}{x_i!},m>0
]
看看上面那条式子，是不是很像？
根据上面那条式子，有：
[[m=1]=sum_{x_1+x_2+...+x_k=m} frac{m!}{k!}prod_{i=1}^kfrac{f(x_i)}{x_i!}\
herefore frac{[m=1]}{m!}=[m=1]=sum_{x_1+x_2+...+x_k=m} frac{1}{k!}prod_{i=1}^kfrac{f(x_i)}{x_i!}\
herefore [x^m]e^{F(x)}=[m=1]
]
另外有特殊的一点，我们让(f(0)=0)，于是就有([x^0]F(x)=0)，所以([x^0]e^{F(x)}=1)，那么(e^{F(x)}=1+x)。
则有：
[F(x)=ln(1+x)\
herefore F'(x)=frac{1}{(1+x)}\
herefore F'(x)=sum_{i=0}^{infty} (-1)^ix^i\
herefore F(x)=sum_{i=0}^{infty} frac{(-1)^{i-1}}{i}x^i\
herefore frac{f(i)}{i!}=frac{(-1)^{i-1}}{i}\
herefore f(i)=(-1)^{i-1}(i-1)!
]
至此，本题推导完毕。
code:
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
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?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int qpow(int n,int k){
int ret=1;
while(k){
if(k&1)ret=mul(ret,n);
n=mul(n,n);
k>>=1;
}
return ret;
}
#define pii pair<int,int>
#define mp make_pair
#define fi first
#define se second
int mxn;
int ss[10010];
int tot;
int prime[10010];
bool vis[10010];
int siz[10010];
int fac[10010],ifac[10010];
int calc(int n,int i){
int ret=0;
while(n){
ret+=n/i;
n/=i;
}
return ret;
}
void Euler(int n){
fac[0]=1;
for(int i=1;i<=n;++i)fac[i]=mul(fac[i-1],i);
ifac[n]=qpow(fac[n],mod-2);
for(int i=n-1;i>=0;--i)ifac[i]=mul(ifac[i+1],i+1);
for(int i=2;i<=n;++i){
if(!vis[i]){
prime[++tot]=i;
ss[i]=calc(n,i);
siz[ss[i]]++;
mxn=max(mxn,ss[i]);
}
for(int j=1;j<=tot&&i*prime[j]<=n;++j){
vis[i*prime[j]]=true;
if(i%prime[j]==0)break;
}
}
}
int n,m;
int sum;
int a[10010];
int h[50][10010];
map<vector<int>,int> divide;
void dfs(int now,int st,int all,vector<int> D){
if(all==sum){
int ret=1;
for(int i=1;i<=mxn;++i){
ret=mul(ret,qpow(h[now][i],siz[i]));
}
divide[D]=ret;
return;
}
for(int i=st;i+all<=sum;++i){
if(i+all!=sum&&2*i+all>sum)continue;
for(int j=0;j<=mxn;++j){
h[now+1][j]=h[now][j];
}
for(int j=i;j<=mxn;++j){
h[now+1][j]=add(h[now+1][j],h[now+1][j-i]);
}
D.push_back(i);
dfs(now+1,i,i+all,D);
D.pop_back();
}
}
map<vector<pii>,int> dp[2];
map<vector<pii>,int>::iterator it;
int ans;
void solve(){
vector<pii> empty;
empty.clear();
int now=0,lst=1;
dp[lst][empty]=1;
for(int i=1;i<=m;++i){
dp[now].clear();
for(it=dp[lst].begin();it!=dp[lst].end();++it){
vector<pii> nw;
for(int j=0,up=(*it).fi.size();j<up;++j){
nw=(*it).fi;
nw[j].fi+=a[i];
nw[j].se++;
sort(nw.begin(),nw.end());
dp[now][nw]=add(dp[now][nw],(*it).se);
}
nw=(*it).fi;
nw.push_back(mp(a[i],1));
sort(nw.begin(),nw.end());
dp[now][nw]=add(dp[now][nw],(*it).se);
}
now^=1,lst^=1;
}
for(it=dp[lst].begin();it!=dp[lst].end();++it){
vector<pii> nw=(*it).first;
vector<int> vec;
vec.clear();
int tmp=((m-nw.size())&1)?dec(mod,1):1;
for(int i=0;i<nw.size();++i){
vec.push_back(nw[i].fi);
tmp=mul(tmp,fac[nw[i].se-1]);
}
ans=add(ans,mul(tmp,mul(divide[vec],(*it).se)));
}
}
int main(){
scanf("%d%d",&n,&m);
Euler(n);
for(int i=1;i<=m;++i){
scanf("%d",&a[i]);
sum+=a[i];
}
h[0][0]=1;
vector<int> empty;empty.clear();
dfs(0,1,0,empty);
sort(a+1,a+1+m);
int sb=1,dv=1;
for(int i=2;i<=m;++i){
if(a[i]==a[i-1])sb++;
else {
dv=mul(dv,ifac[sb]);
sb=1;
}
}
dv=mul(dv,ifac[sb]);
solve();
printf("%d
",mul(ans,dv));
}