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

那么，对于每一种方案，它的贡献就是(分配质因数的方案数 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));
}