【XSY3657】 因数分解

【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来解决。

那么，对于每一种方案，它的贡献就是\(分配质因数的方案数\times 划分a_i的方案数\times 容斥系数\)。

于是来考虑容斥系数。

不妨大胆假设，容斥系数只和划分到同一组的\(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\)，则有：

\[\because F(x)=\sum_{i=0}^{\infty} \frac{f(i)}{i!}x^i\\ \therefore 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!}\\ \therefore [x^m]e^{F(x)}=\sum_{x_1+x_2+...+x_k=m} \frac{1}{k!}\prod_{i=1}^k\frac{f(x_i)}{x_i!},m>0 \]

看看上面那条式子，是不是很像？

根据上面那条式子，有：

\[[m=1]=\sum_{x_1+x_2+...+x_k=m} \frac{m!}{k!}\prod_{i=1}^k\frac{f(x_i)}{x_i!}\\ \therefore \frac{[m=1]}{m!}=[m=1]=\sum_{x_1+x_2+...+x_k=m} \frac{1}{k!}\prod_{i=1}^k\frac{f(x_i)}{x_i!}\\ \therefore [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)\\ \therefore F'(x)=\frac{1}{(1+x)}\\ \therefore F'(x)=\sum_{i=0}^{\infty} (-1)^ix^i\\ \therefore F(x)=\sum_{i=0}^{\infty} \frac{(-1)^{i-1}}{i}x^i\\ \therefore \frac{f(i)}{i!}=\frac{(-1)^{i-1}}{i}\\ \therefore 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\n",mul(ans,dv));
}