考慮這題是什麼意思，其實就是讓你把 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;
}