(color{#0066ff}{题解 })
可以发现， 数据范围中的n特别小，容易想到状压
可以想到类似于状压DP的思路，按列进行转移
那么应该有3维，(f[i][j][k])代表到第i列，j的每一位表示这一行有多少连续的男生，k表示当前有多少列全是男生，的方案数
看到m的范围，我们肯定是要找一个(O(logm))的东西加速转移，自然是矩阵加速
然后我们来看看有多少个状态，看看是否可行
j有(p^n)个，k有q个（用矩阵转移第一维自然不需要）
那么状态数依然达到了一个(p^n*q=3^8*3=19683)
这东西要是(O(n^3))去转移绝对T到飞起，而且空间也是开不下
所以我们考虑减少状态的数量
发现，我们并不需要知道每一行当前有多少个连续的男生，我们只需要知道有多少行的连续的男生是(0,1dots p-1)个
这相当于把一个数n拆分成p个数的和
这样的状态数绝对是少得多了，就可以矩阵加速了
我们可以通过dfs预处理出所有的状态（状态数极限$le$200）
然后为了推出转移矩阵，我们需要考虑一个状态对另一个状态的贡献，直接(n^2)枚举即可
假设我们现在考虑(x o y)的贡献
设(x_0, x_1, x_2,xx)表示状态x中，有连续0个男生的行的数量，有连续1个男生的行的数量，有连续2个男生的行的数量，有多少列全是男生，这个已经dfs预处理过，y同理
如果(xx + 1 == yy)
那就说明当前这一列必须全是男生，于是(x_0)要全给(y_1)，(x_1)要全给(y_2)，而且(x_2)必须为0，这样才有1的贡献
然后考虑不全为1的
我们肯定是放女生或男生，那么拿(x_1)来说，如果一些行放女生，那么这些行就要转移到(y_0)，如果放男生，就要转移到(y_2)，那么显然组合一下，方案就是(C_{x_1}^{x_1-y_2})，(x_0)同理
然后构造完转移矩阵，加速递推即可
#include<bits/stdc++.h>
#define LL long long
LL in() {
char ch; LL x = 0, f = 1;
while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
return x * f;
}
const int mod = 1e9 + 7;
LL n, m, p, q, len;
struct Matrix {
LL ju[222][222];
Matrix() { memset(ju, 0, sizeof ju); }
void e() {
memset(ju, 0, sizeof ju);
for(int i = 1; i <= len; i++) ju[i][i] = 1;
}
friend Matrix operator * (const Matrix &a, const Matrix &b) {
Matrix c;
for(int i = 1; i <= len; i++)
for(int j = 1; j <= len; j++)
for(int k = 1; k <= len; k++)
(c.ju[i][j] += a.ju[i][k] * b.ju[k][j] % mod) %= mod;
return c;
}
Matrix ksm(LL y) {
Matrix x = *this, re;
re.e();
while(y) {
if(y & 1) re = re * x;
x = x * x;
y >>= 1;
}
return re;
}
}beg, A;
LL fac[22], inv[22];
int ls[22];
struct node {
int st[3], num;
node(int a = 0, int b = 0, int c = 0, int d = 0): num(d) {
st[0] = a, st[1] = b, st[2] = c;
}
}st[105050];
LL ksm(LL x, LL y) {
LL re = 1LL;
while(y) {
if(y & 1) re = re * x % mod;
x = x * x % mod;
y >>= 1;
}
return re;
}
LL C(int x, int y) {
return fac[x] * inv[y] % mod * inv[x - y] % mod;
}
void dfs(int dep, int tot) {
if(dep == p + 1) {
if(tot) return;
for(int i = 0; i <= q; i++) st[++len] = node(ls[1], ls[2], ls[3], i);
return;
}
for(int i = tot; i >= 0; i--) ls[dep] = i, dfs(dep + 1, tot - i);
}
void predoit() {
fac[0] = inv[0] = 1;
for(int i = 1; i <= 20; i++) fac[i] = 1LL * i * fac[i - 1] % mod, inv[i] = ksm(fac[i], mod - 2);
}
void work(const node &a, const node &b, int x, int y) {
int a1 = a.st[0], b1 = a.st[1], c1 = a.st[2], d1 = a.num;
int a2 = b.st[0], b2 = b.st[1], c2 = b.st[2], d2 = b.num;
if(d1 + 1 == d2) {
if(a1 == b2 && b1 == c2 && !c1) return (void)(A.ju[y][x]++);
} else if(d1 == d2) {
if(a1 == b2 && b1 == c2 && !c1) return;
if(a1 >= b2 && b1 >= c2) return (void)(A.ju[y][x] = C(a1, a1 - b2) * C(b1, b1 - c2) % mod);
}
}
void Build_Matrix() {
dfs(1, n);
#ifdef olinr
for(int i = 1; i <= len; i++) printf("%d %d %d %d
", st[i].st[0], st[i].st[1], st[i].st[2], st[i].num);
#endif
for(int i = 1; i <= len; i++)
for(int j = 1; j <= len; j++)
work(st[i], st[j], i, j);
beg.ju[1][1] = 1;
}
void query() {
beg = A.ksm(m) * beg;
LL ans = 0;
for(int i = 1; i <= len; i++) (ans += beg.ju[i][1]) %= mod;
printf("%lld", ans);
}
int main() {
freopen("b.in", "r", stdin);
freopen("b.out", "w", stdout);
n = in(), m = in(), p = in(), q = in();
predoit();
Build_Matrix();
query();
return 0;
}