題目大意

給定n,k，求把n!拆分成k個不同的正整數的乘積的方案數。（一種方案的排列仍是一種方案）。答案對109+7取模。

n≤10000 k≤30

時限為4s

分析

這是一道容斥好題。

首先可以不管算重，最終答案除以k!即可。

接下來考慮如何容斥。k個數互不相同，其實就相當於k(k−1)2個限制條件。如果一個方案不滿足x個條件，那麼它要乘上的容斥係數為(−1)x

如果直接列舉不滿足哪些，顯然會超時。不妨轉換一下：已知條件是形如(x,y)的，表示第x個數≠第y個數。把k個數看成k個點，限制(x,y)看成連線x,y的邊。那麼一個選擇方案相當於選擇一些邊，使其變成許多個聯通塊。

好了，現在可以嘗試列舉每個聯通塊的大小。這就相當於對k的整數拆分。為了避免算重，當前列舉的這一個數必須不小於前一個。跑出來，30能有5000多種拆法。那麼整數拆分的方案數就可以接受了。接下來到下一步。

現在已經確定了k的一個整數拆分方案（即若干個聯通塊大小）。假設有m個聯通塊，第i個大小為si。它對答案的貢獻是多少？這要考慮三個東西：
1. 把k個點分配到這些聯通塊的方案數。
2. 對應的分配邊的所有方案對答案的貢獻
3. 拆分n!，滿足形成這些聯通塊的方案數。
它對答案的貢獻就是三個數的乘積。現在一個個來考慮。

首先是第1部分

這個比較簡單。首先總的排列數有k!個。設在當前的整數拆分方案中，數字i出現了pi次，即∑i∗pi=k。那麼對於拆出來的每個聯通塊si，這表示有si個數相同，就要除掉si!，接下來由於i出現了pi次，又要除掉pi!。最終得到：

第2部分

對於每個聯通塊分別考慮。因為一條邊的意義是兩個數相同，所以一個聯通塊裡的所有數都是相等的。那就可以直接考慮容斥係數的總和了。
設S(n)表示n個點的聯通圖的集合，e(G)表示圖G的邊數。對於第i個聯通塊（大小為si），它對答案的貢獻就是∑G∈S(si)(−1)e(G)。
設g(n)等於n個點的聯通圖對答案的貢獻。這個可以考慮用容斥計算，即全集減去不合法方案數。
首先考慮不合法的。其實就是n個點的不聯通圖。它形成了至少2個聯通塊。列舉編號為1的點所在的聯通塊大小，那麼：

考慮f(n)的值。顯然f(1)=1。當n>1，就會有n(n−1)2條可能出現的邊。如果其中一條邊出現了，容斥係數就乘-1，否則乘1。所有方案加起來後發現，當n>1時，f(n)=0！
上面的式子變成：

最後是第3部分

首先給n!分解質因數。每個質因子都是相對獨立的，就可以分開考慮。
對於一個質因子，設它共出現了cnt次，現在要把它分配到m個聯通塊裡。一個聯通塊裡的數要滿足相同，每個質因子出現個數也一定要相同。所以對於si，它要在cnt份裡佔掉si的倍數份（可以是0）。這就是個完全揹包問題了。

優化

為了提高執行效率，可能需要一些優化。
1. 遞迴對k進行拆分，然後做揹包時相當於有m個物品。回溯的時候，前m-1個物品的答案要保留，不用每個方案都重新做一次dp
2. 假設k當前剩下x沒有拆分，先列舉到x/2，再直接整個拆分
3. 適當去掉一些模運算

我的程式極限資料跑了3s左右

#include <cstdio>
#include <cstring>
#include <algorithm>

#define Min(a,b) ((a)<(b)?(a):(b))

using namespace std;

const int N=10005,M=31,mo=1e9+7,NN=100000;

typedef long long LL;

int n,m,Fac[N],Inv[N],F_Inv[N],g[N],tot,p[N],f[M][NN],cnt[N],st[N],sum,ans;

bool bz[N];

void init()
{
    scanf("%d%d",&n,&m);
    for (int i=2;i<=n;i++)
    {
        if (!bz[i]) p[tot++]=i;
        for (int j=0;j<tot && i*p[j]<=n;j++)
        {
            bz[i*p[j]]=1;
            if (i%p[j]==0) break;
        }
    }
    for (int i=2;i<=n;i++)
        for (int j=0,x=i;j<tot;j++)
            for (;x%p[j]==0;x/=p[j]) cnt[j]++;
    for (int i=1;i<=tot;i++) st[i]=st[i-1]+cnt[i-1]+1;
    Fac[0]=F_Inv[0]=Fac[1]=Inv[1]=F_Inv[1]=1;
    for (int i=2;i<=n;i++)
    {
        Fac[i]=(LL)Fac[i-1]*i%mo;
        Inv[i]=(LL)Inv[mo%i]*(mo-mo/i)%mo;
        F_Inv[i]=(LL)F_Inv[i-1]*Inv[i]%mo;
    }
}

void dfs(int x,int y,int m,int sum,int t)
{
    if (m==0)
    {
        for (int i=0;i<tot;i++) sum=(LL)sum*f[x-1][st[i+1]-1]%mo;
        ans=(ans+sum)%mo;
        return;
    }
    if (y>m) return;
    for (int i=0;i<tot;i++)
    {
        for (int j=st[i];j<st[i]+y;j++) f[x][j]=f[x-1][j];
        for (int j=st[i]+y;j<st[i+1];j++) f[x][j]=(f[x-1][j]+f[x][j-y])%mo;
    }
    dfs(x+1,y,m-y,(LL)sum*F_Inv[y]%mo*Inv[t+1]%mo*g[y]%mo,t+1);
    if (y==m) return;
    for (int k=y+1;k*2<=m;k++)
    {
        for (int i=0;i<tot;i++)
        {
            f[x][st[i]+k-1]=f[x-1][st[i]+k-1];
            for (int j=st[i]+k;j<st[i+1];j++) f[x][j]=(f[x-1][j]+f[x][j-k])%mo;
        }
        dfs(x+1,k,m-k,(LL)sum*F_Inv[k]%mo*g[k]%mo,1);
    }
    for (int i=0;i<tot;i++)
    {
        for (int j=st[i];j<st[i]+m;j++) f[x][j]=f[x-1][j];
        for (int j=st[i]+m;j<st[i+1];j++) f[x][j]=(f[x-1][j]+f[x][j-m])%mo;
    }
    dfs(x+1,m,0,(LL)sum*F_Inv[m]%mo*g[m]%mo,1);
}

void work()
{
    g[1]=1;
    for (int i=2;i<=m;i++) g[i]=(-(LL)(i-1)*g[i-1])%mo;
    ans=0;
    for (int i=0;i<tot;i++) f[0][st[i]]=1;
    dfs(1,1,m,Fac[m],0);
    ans=(LL)(ans+mo)%mo*F_Inv[m]%mo;
    printf("%d\n",ans);
}

int main()
{
    //freopen("jdt.in","r",stdin); freopen("jdt.out","w",stdout);
    init();
    work();
    fclose(stdin); fclose(stdout);
    return 0;
}