【GDOI2020模拟01.17】小 ω 玩游戏 （容斥+EGF）

小 ω 正在玩一个游戏。
小 ω 有一个 n 行 m 列的网格，初始每个方格中都有数字 0。她需要执行 q 次操作，每次操作可以选择其中一个方格 (x, y)，然后先将第 x 行的数全部 +1，接着将第 y 列的数全部 +1。
小 ω 想知道有多少种执行操作的方式能使最后的网格中有不超过 k 个奇数。
两种方式不同当且仅当存在某一步中选择的方格坐标不同。

(1<=n,m<=2e5,q<=10^{18})

---

考虑行列分开，对行，算出(f(x))表示恰好x行奇数的方案数，对列同理，算出设为g。

有了f和g就可以直接解个不等式然后前缀和统计。

考虑有n行，有p行奇数的方案数，直接列出EGF：

(f[p]=({e^x-e^{-x}over 2})^p*({e^x+e^{-x}over 2})^{n-p}[x^q]*q!*C_n^p)

暴力展开是(O(n^3))的。

考虑先求出(({e^x-e^{-x}over 2})^n)，然后每次除除乘乘就可以做到(O(n^2))的。

然后优化不了了，此时应该考虑容斥：

至少p个奇数行的EGF就很简单了：

(f2(p)=({e^x-e^{-x}over 2})^p*e^{x(n-p)}[x^q]*q!*C_n^p)

(=sum_{i=0}^p C_{p}^i*(-1)^{p-i}*e^{xi-x(p-i)}/2^p*e^{x(n-p)}[x^q]*q!*C_n^p)

(=sum_{i=0}^p C_p^i*(-1)^{p-i}*e^{(n-2(p-i))x}[x^q]/*q!/2^p*C_n^p)

(=sum_{i=0}^pC_p^i*(-1)^{p-i}*(n-2(p-i))^q/2^p*C_{n}^p)

显然可以卷积求。

然后有(f2[n]=sum_{i>=n}f[i]*C_{i}^n)

所以(f[n]=sum_{i>=n}f2[i]*C_{i}^n*(-1)^{i-n})

再做一个卷积就可以求出f。

#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;

const int mo = 998244353;

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;
}

const int nm = 1 << 19;

int r[nm]; ll w[nm];
ll a[nm], b[nm];

void build() {
	for(int i = 1; i < nm; i *= 2) {
		w[i] = 1; ll v = ksm(3, (mo - 1) / 2 / i);
		ff(j, 1, i) w[i + j] = w[i + j - 1] * v % mo;
	}
}
void dft(ll *a, int n, int f) {
	ff(i, 0, n) {
		r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
		if(i < r[i]) swap(a[i], a[r[i]]);
	} ll b;
	for(int i = 1; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i) ff(k, 0, i)
		b = a[i + j + k] * w[i + k], a[i + j + k] = (a[j + k] - b) % mo, a[j + k] = (a[j + k] + b) % mo;
	if(f == -1) {
		reverse(a + 1, a + n);
		b = ksm(n, mo - 2);
		ff(i, 0, n) a[i] = (a[i] + mo) * b % mo;
	}
}
void fft(ll *a, ll *b, int n) {
	dft(a, n, 1); dft(b, n, 1);
	ff(i, 0, n) a[i] = a[i] * b[i] % mo;
	dft(a, n, -1);
}

const int N = 2e5 + 5;

int n, m; ll q, k;
ll f[N], g[N], s1[N], s2[N];

ll fac[N], nf[N];

ll s[N];

ll C(int n, int m) {
	return fac[n] * nf[m] % mo * nf[n - m] % mo;
}

void calc(ll *f, int n) {
	fac[0] = 1; fo(i, 1, n) fac[i] = fac[i - 1] * i % mo;
	nf[n] = ksm(fac[n], mo - 2); fd(i, n, 1) nf[i - 1] = nf[i] * i % mo;
	
	int tt = 0;
	while((1 << (++ tt)) <= 2 * n);
	
	memset(a, 0, sizeof a);
	memset(b, 0, sizeof b);
	fo(i, 0, n) a[i] = nf[i], b[i] = (i & 1 ? -1 : 1) * ksm(n - 2 * i, q) * nf[i] % mo;
	fft(a, b, 1 << tt);
	ll ni2 = ksm(2, mo - 2), s = 1;
	fo(i, 0, n) a[i] = a[i] * fac[i] % mo * s % mo * C(n, i) % mo, s = s * ni2 % mo;
	ff(i, n + 1, 1 << tt) a[i] = 0;
	memset(b, 0, sizeof b);
	fo(i, 0, n) a[i] = a[i] * fac[i] % mo, b[n - i] = nf[i] * (i & 1 ? -1 : 1);
	fft(a, b, 1 << tt);
	fo(i, 0, n) f[i] = a[i + n] * nf[i] % mo;
}

void End() {
	s2[0] = g[0]; fo(i, 1, m) s2[i] = (s2[i - 1] + g[i]) % mo;
	ll ans = 0;
	for(ll x = 0; x <= n; x ++) {
		int xs = n - 2 * x;
		if(xs == 0) {
			if(x * m <= k) ans += f[x] * s2[m] % mo;
		} else
		if(xs > 0) {
			ll t = floor((long double) (k - x * m) / xs);
			t = min(t, (ll) m);
			if(t >= 0) ans += f[x] * s2[t] % mo;
		} else {
			ll t = ceil((long double) (k - x * m) / xs);
			if(t <= m) ans += f[x] * (t > 0 ? (s2[m] - s2[t - 1]) : s2[m]) % mo;
		}
	}
	ans = (ans % mo + mo) % mo;
	pp("%lld
", ans);
}

int main() {
	freopen("play.in", "r", stdin);
	freopen("play.out", "w", stdout);
	build();
	scanf("%d %d %lld %lld", &n, &m, &q, &k);
	calc(f, n); calc(g, m);
	End();

}