考虑这题是什么意思，其实就是让你把 DAG 划分成若干个集合，点之间连边转化为对应集合之间连边以后图仍然是一个 DAG，然后需要知道划分成了多少个集合，每种集合的个数求出方案数，乘上对应的系数并求和。
系数是很显然的，即：
\[{k+1\choose i}\frac{i!k!}{n!\prod_{i=1}^k (n+i)} \]考虑怎么求方案数。记 \(f_{S,i}\) 表示把集合 \(S\) 分成 \(i\) 个子集且最终的图是一个 DAG 的方案数。考虑枚举一个没有出边的子集然后转移，这样可以保证最终的图是一个 DAG，但是会算重。
会算重，所以就可以做一个容斥。记 \(g_{S,i}\) 表示把集合 \(S\) 分成 \(i\) 个子集且这些子集之间无边的方案数。可以直接枚举子集然后转移，使得子集中包含 \(S\) 中编号最小的点，这样求出来的方案数是正确的。
然后就是最终的容斥部分，易得：
\[f_{S,i}=\sum_{S'\subseteq S}\sum_{j=1}^{i}(-1)^{j-1}\cdot g_{S',j}\cdot f_{S/S',i-j} \]最后乘上系数求和即可。时间复杂度
\(T(n)=\sum_{i=0}^n{n\choose i}2^ii^2\mathcal{O}(1)=2\times 3^{n-2}n(2n+1)\mathcal{O}(1)=\mathcal{O}(3^nn^2)\)，常数很小，可以通过此题。
参考代码：
#include<bits/stdc++.h>
#define ll long long
#define mxn 200003
#define md 1000000007
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define rept(i,a,b) for(int i=a;i<b;++i)
using namespace std;
inline ll power(ll x,int y){
ll ans=1;
for(;y;y>>=1){
if(y&1)ans=ans*x%md;
x=x*x%md;
}
return ans;
}
int n,m,k,t,d[18],d1[18],dd[1<<15],s[18],a[18],c[20][20];
ll xi,ans,fac[20],dp[1<<15],f1[1<<15][17],f[1<<15][17];
signed main(){
scanf("%d%d%d",&n,&m,&k);
c[0][0]=1;
rep(i,1,18){
c[i][0]=1;
rep(j,1,i)c[i][j]=(c[i-1][j-1]+c[i-1][j])%md;
}
fac[0]=1;
rep(i,1,18)fac[i]=fac[i-1]*i%md;
for(int i=0,x,y;i<m;++i){
scanf("%d%d",&x,&y);
x--,y--;
d[y]|=1<<x,d1[x]|=1<<y;
}
xi=1;
rep(i,1,n)xi=xi*i%md;
xi=power(xi,md-2);
rep(i,1,k)xi=xi*power(n+i,md-2)%md*i%md;
dp[0]=1;
rept(s,1,1<<n){
rept(i,0,n)if((s>>i)&1){
if(d[i]&s)continue;
dp[s]=(dp[s]+dp[s^(1<<i)])%md;
}
}
rept(s,1,1<<n){
rept(i,0,n)if((s>>i)&1){
dd[s]=dd[s^(1<<i)]|d[i];
break;
}
}
f1[0][0]=1;
rept(s,1,1<<n){
int s1=s^(s&-s),ct=min(__builtin_popcount(s),k+1);
for(int s2=s1;;s2=(s2-1)&s1){
rept(i,0,n)if(((s2|(s&-s))>>i)&1&&((d[i]|d1[i])&(s1^s2)))goto nxt;
rep(i,1,ct)f1[s][i]=(f1[s][i]+f1[s1^s2][i-1]*dp[s2|(s&-s)])%md;
nxt:;
if(!s2)break;
}
}
f[0][0]=1;
rept(s,1,1<<n){
int ct=min(__builtin_popcount(s),k+1);
for(int s1=s;s1;s1=(s1-1)&s){
if((s^s1)&dd[s1])continue;
rep(i,1,ct){
rep(j,1,i){
if(j&1)f[s][i]=(f[s][i]+f[s^s1][i-j]*f1[s1][j])%md;
else f[s][i]=(f[s][i]-f[s^s1][i-j]*f1[s1][j])%md;
}
}
}
}
rep(i,1,k+1)ans=(ans+f[(1<<n)-1][i]*c[k+1][i]%md*fac[i])%md;
cout<<(ans*xi%md+md)%md;
return 0;
}