分析:
又被开除人籍了。。
首先我们要把行和列分开计算
假设\(n\)行中有\(i\)个被访问了奇数次,\(m\)列中有\(j\)个被访问了奇数次
那么最终的奇数点数可以计算为\((n-i)j+(m-j)i\)
(然而这步我都没想到,7分暴力走人。。。
于是我们把行列分开处理做\(q\)次
我们现在只看行,列做同样处理
设\(f_i\)表示\(n\)行做了\(q\)次里面恰好有\(i\)行为奇数
看样子像是要用生成函数去表示它
啥啊,我是废物不会啊
考虑容斥,设\(g_i\)表示\(n\)行做了\(q\)次里面至少有\(i\)行为奇数
关系式为:
\(g_i=\sum_{j=i}^{n}f_i\)
二项式反演:
\(f_i=\sum_{j=i}^{n}(-1)^{j-i}g_i\)
我们知道了\(g\),这个卷积形式喜闻乐见NTT就可以处理出\(f\)
来算\(g\)
先把\(i\)个必为奇数的放到前面,最后系数乘一个\(\binom{n}{i}\)就好
表示奇数位并带阶乘逆元的生成函数是\(\frac{e^x-e^{-x}}{2}\)
所以
\[g_i=\binom{n}{i}q![x^q](\frac{e^x-e^{-x}}{2})^i(e^x)^{n-i}
\]
取第\(q\)项,前面的\(q!\)会和生成函数第\(q\)项的阶乘逆元抵消
用一下二项式定理:
\[g_i=\binom{n}{i}2^{-i}\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}q![x^q]e^{(n−2i+2j)x}
\]
\[g_i=\binom{n}{i}2^{-i}\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}(n−2i+2j)^q
\]
喜闻乐见NTT可以求\(g\)
之后再用NTT求\(f\)
行和列都得到了各自的\(f\)
我们设一个为\(f\)一个为\(g\)吧
考虑式子:
\[Ans=\sum_{(n-i)j+(m-j)i\leq k}f_ig_j
\]
我们枚举\(i\),合法的\(j\)是一个区间
那么我们对\(g\)求一个前缀和,分类讨论计算即可
讨论起来有点恶心,而且注意开long long
复杂度\(O(nlogn)\)
#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<set>
#include<map>
#include<vector>
#include<string>
#define maxn 1000005
#define INF 0x3f3f3f3f
#define MOD 998244353
using namespace std;
inline long long getint()
{
long long num=0,flag=1;char c;
while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;
while(c>='0'&&c<='9')num=num*10+c-48,c=getchar();
return num*flag;
}
int n,m,len;
long long q,k;
int a[maxn],b[maxn];
int f[maxn],g[maxn],sum[maxn];
int rev[maxn];
int fac[maxn],inv[maxn];
inline int upd(int x){return x<MOD?x:x-MOD;}
inline int C(int p,int q)
{return 1ll*fac[p]*inv[q]%MOD*inv[p-q]%MOD;}
inline int ksm(int num,int k)
{
int ret=1;
for(;k;k>>=1,num=1ll*num*num%MOD)if(k&1)ret=1ll*ret*num%MOD;
return ret;
}
inline void NTT(int *a,int N,int opt)
{
for(int i=0;i<N;i++)if(rev[i]<i)swap(a[i],a[rev[i]]);
for(int i=1;i<N;i<<=1)
{
int wn=ksm(3,(MOD-1)/(i<<1));
if(!~opt)wn=ksm(wn,MOD-2);
for(int j=0;j<N;j+=i<<1)for(int k=0,w=1;k<i;k++,w=1ll*w*wn%MOD)
{
int x=a[j+k],y=1ll*w*a[i+j+k]%MOD;
a[j+k]=upd(x+y),a[i+j+k]=upd(x-y+MOD);
}
}
if(!~opt)for(int i=0,Inv=ksm(N,MOD-2);i<N;i++)a[i]=1ll*a[i]*Inv%MOD;
}
inline void solve(int n,int *f)
{
for(int i=0;i<len;i++)a[i]=b[i]=0;
for(int i=0;i<=n;i++)a[i]=inv[i];
for(int i=0;i<=n;i++)b[i]=1ll*(i&1?MOD-1:1)*inv[i]%MOD*ksm(upd(n-2*i+MOD),q%(MOD-1))%MOD;
NTT(a,len,1),NTT(b,len,1);
for(int i=0;i<len;i++)a[i]=1ll*a[i]*b[i]%MOD;
NTT(a,len,-1);
for(int i=0;i<=n;i++)f[i]=1ll*a[i]*ksm(inv[2],i)%MOD*fac[n]%MOD*inv[n-i]%MOD;
for(int i=0;i<len;i++)a[i]=b[i]=0;
for(int i=0;i<=n;i++)a[i]=1ll*f[i]*fac[i]%MOD;
for(int i=0;i<=n;i++)b[n-i]=1ll*(i&1?MOD-1:1)*inv[i]%MOD;
NTT(a,len,1),NTT(b,len,1);
for(int i=0;i<len;i++)a[i]=1ll*a[i]*b[i]%MOD;
NTT(a,len,-1);
for(int i=0;i<=n;i++)f[i]=1ll*a[i+n]*inv[i]%MOD;
}
int main()
{
n=getint(),m=getint(),q=getint(),k=getint();
if(!n||!m){printf("0\n");return 0;}
if(n<m)swap(n,m);
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=2;i<=n;i++)fac[i]=1ll*fac[i-1]*i%MOD;
for(int i=2;i<=n;i++)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
for(int i=2;i<=n;i++)inv[i]=1ll*inv[i]*inv[i-1]%MOD;
len=1;
while(len<=2*n)len<<=1;
for(int i=0;i<len;i++)rev[i]=(rev[i>>1]>>1)|(i&1?len>>1:0);
solve(n,f),solve(m,g);
for(int i=1;i<=m;i++)g[i]=upd(g[i-1]+g[i]);
int ans=0;
for(int i=0;i<=n;i++)
{
if(n>2*i)
{
long long p=floor((k-1.0*i*m)/(n-2*i));
if(p>=0)ans=upd(ans+1ll*f[i]*g[min(1ll*m,p)]%MOD);
}
else if(n<2*i)
{
long long p=ceil((k-1.0*i*m)/(n-2*i));
if(p<=m)ans=upd(ans+1ll*f[i]*upd(g[m]-(p>0?g[p-1]:0)+MOD)%MOD);
}
else if(k-1ll*i*m>=0)ans=upd(ans+1ll*f[i]*g[m]%MOD);
}
printf("%d\n",ans);
}