题目大意：
(n)种花色，每个花色(m)种数字，每个花色数字牌不超过(4)张的雀魂。
随机选(3k+2)个牌，问能拆成(k)个面子和(1)个对子的概率。
(n,m le 10^9, k le 30)
题解：
先考虑一种花色的牌选出(0..3k+2)的方案数，后面做一个背包就可以组合。
再考虑一种花色拆成若干段(>0)的连续段会比较好做，求出每个连续段（每个数字牌数都(>0)）的方案数，再用一个(O(k^5))的dp也能组合。
问题是求(f[i][j])表示长度(=i)的连续段，用了(j)张牌且合法的方案数。
考虑怎么判断一个段是否合法，可以dp，([i=0 sim 4][j = 0 sim 4][k = 0 sim 1])表示最后两个的剩余大小和有没有选对子。
计数不难想到dp套dp，但是(2^{5*5*2})有点大，可是经过搜索发现有用状态只有(248)个。
于是就做完了，挺好写的。
比赛时忘记搜索缩状态了，GG
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i < _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
ll m, n, k;
const int N = 100;
const int mo = 1e9 + 7;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
ll cnt[N][N];
ll c[5][5];
void add(ll &x, ll y) {
(x += y) %= mo;
}
namespace sub1 {
ll a2[51];
int id[5][5][2], id0;
map<ll, int> bz; int bz0;
int bx[51];
const int M = 300;
int to[M][5];
int ok1[M], ok2[M];
int dg(ll s) {
if(bz[s]) return bz[s];
bz[s] = ++ bz0;
ok1[bz0] = s >> id[0][0][0] & 1;
ok2[bz0] = s >> id[0][0][1] & 1;
fo(x, 1, 4) {
fo(i, 0, 49) bx[i] = 0;
fo(i, 0, 4) fo(j, 0, 4) fo(k, 0, 1) if(s >> id[i][j][k] & 1) {
fo(nk, 0, !k) fo(c, 0, 1) {
int a[4];
a[1] = i, a[2] = j, a[3] = x;
while(a[1] && a[2] && a[3]) a[1] --, a[2] --, a[3] --;
if(a[1]) continue;
a[3] -= nk * 2 + c * 3;
if(a[3] < 0) continue;
bx[id[a[2]][a[3]][k | nk]] = 1;
}
}
ll ns = 0;
fo(i, 0, 49) if(bx[i]) ns |= a2[i];
to[bz[s]][x] = dg(ns);
}
return bz[s];
}
ll f[N][N][M];
void work() {
fo(i, 0, 4) {
c[i][0] = 1;
fo(j, 1, i) c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mo;
}
a2[0] = 1; fo(i, 1, 50) a2[i] = a2[i - 1] * 2;
fo(i, 0, 4) fo(j, 0, 4) fo(k, 0, 1)
id[i][j][k] = ++ id0;
dg(a2[id[0][0][0]]);
f[0][0][1] = 1;
fo(i, 0, 3 * k + 2) {
fo(j, 0, 3 * k + 2) {
fo(u, 1, bz0) if(f[i][j][u]) {
if(ok1[u]) add(cnt[i][j], f[i][j][u]);
if(ok2[u]) add(cnt[i][j], f[i][j][u]);
fo(x, 1, 4) if(j + x <= 3 * k + 2) {
add(f[i + 1][j + x][to[u][x]], f[i][j][u] * c[4][x]);
}
}
}
}
}
}
ll g[N];
ll inv[N];
ll C(int n, int m) {
if(n < m) return 0;
ll s = 1;
fo(i, n - m + 1, n) s = s * i % mo;
fo(i, 1, m) s = s * inv[i] % mo;
return s;
}
namespace sub2 {
ll f[N][N][N];
void work() {
fo(i, 1, 3 * k + 2) inv[i] = ksm(i, mo - 2);
f[0][0][0] = 1;
fo(i, 0, 3 * k + 2) {
fo(j, 0, 3 * k + 2) {
fo(u, 0, 3 * k + 2) if(f[i][j][u]) {
add(g[u], f[i][j][u] * C(m - j + 1, i));
fo(p, 1, 3 * k + 2) fo(q, 1, 3 * k + 2) if(cnt[p][q] && j + p <= 3 * k + 2 && u + q <= 3 * k + 2) {
if(u % 3 == 2 && q % 3 == 2) continue;
add(f[i + 1][j + p][u + q], f[i][j][u] * cnt[p][q]);
}
}
}
}
}
}
namespace sub3 {
ll f[N][N];
void work() {
f[0][0] = 1;
ll ans = 0;
fo(i, 0, 3 * k + 2) {
fo(j, 0, 3 * k + 2) if(f[i][j]) {
fo(u, 1, 3 * k + 2) if(j + u <= 3 * k + 2 && g[u]) {
if(j % 3 == 2 && u % 3 == 2) continue;
add(f[i + 1][j + u], f[i][j] * g[u]);
}
}
ans = (ans + C(n, i) * f[i][3 * k + 2]) % mo;
}
pp("%lld
", ans);
}
}
int main() {
freopen("majsoul.in", "r", stdin);
freopen("majsoul.out", "w", stdout);
scanf("%lld %lld %lld", &m, &n, &k);
sub1 :: work();
sub2 :: work();
sub3 :: work();
}