Codechef Oct chanllenge Queries on Matrix-JIIT

首先发现矩阵的两个维度显然是互不相干的，假设最后操作后有$x$列被操作奇数次，$y$行操作奇数次

那么最后为奇数的格子个数就是$x(m-y)+(n-x)y$

考虑求出$q$操作后有$x$个位置被操作奇数次的方案数

考虑一个Naive的$dp$，令$dp_{i,j}$为操作$i$次后有$j$位置操作奇数的方案数，显然得到转移为

$dp_{i,j}\cdot j\rightarrow dp_{i+1,j-1}$

$dp_{i,j}\cdot (n-j)\rightarrow dp_{i+1,j+1}$

直接$dp$复杂度为$O(nq)$，用矩阵优化复杂度为$O(n^3\log q)$

没有考虑过求向量列的现行递推式？说不定暴力求递推然后。。。

$O(n^3)$

由于笔者不会数学，所以考虑一个非常暴力非常直观的理解，可以完全抛开组合意义

每次是挑选一个位置异或上1，用形式幂级数可能比较蛋疼，不如直接搞成集合幂级数

即令集合$S$包含所有被操作奇数次的位置，用一个多项式$F=\sum_S a_S x^S$表示答案

那么转移多项式即为$F=\sum_{i=1}^{n} x^{\lbrace i \rbrace}$，转移运算为集合对称差运算(哎就是异或)

那么实际就是要求出$F^q$，可以直接用$\text{FWT}$优化，先$\text{FWT}$，然后求出每一项的$q$次幂，然后$\text{FWT}$回来

(这样不是$n2^n$的吗)

显然的可以发现，$F^i$的任意位置系数$[x^S]F^i$只与$|S|$有关，所以实际上只需要存下含有$0,1,2\cdots,n$个元素的项的系数即可

考虑在这样的多项式上模拟原先的$\text{FWT}$过程

按照快速沃尔什变换的式子，令$G=\text{FWT(F)}=\sum_{S}x^S \sum_{T}(-1)^{|S\cap T|}\cdot [x^T]F $

考虑枚举$|S|,|T|,|S\cap T|$，然后组合数计算系数，令$F_i$表示$[x^S]F(|S|=i)$

则有$\begin{aligned} G_i=\sum_{j}F_j\sum_k C_i^k\cdot C_{n-i}^{j-k}(-1)^{k}\end{aligned}$，其中$C_{n-i}^{j-k}$表示从不相交的部分里选出$j$中剩下的元素

按照该式即可完成$O(n^3)$模拟$\text{FWT}$，注意$\text{IFWT}$时需要除掉$2^n$

算上快速幂复杂度应为$O(n^3+n\log q)$

$O(n^2\log n)$

用$\text{NTT}$优化上式

$O(n^2)$

依然考虑上面$\text{FWT}$的转移式，发现$i$项得到$j$项的贡献为一个常数，考虑直接计算这个常数$W_{i,j}=\sum_k C_i^k\cdot C_{n-i}^{j-k}(-1)^{k}$

我们知道组合数就是二项展开的结果，所以发现实际就是$W_{i,j}=[x^j] (-x+1)^i\cdot (x+1)^{n-i}$

该式可以$O(n^2)$按照$i$递推求出，每次乘上一个$\frac{1-x} {1+x}$即可

实际复杂度为$O(n^2+n\log q)$

实际上应该还可以处理一些稍微复杂点的问题，比如每次可以操作若干个位置

不知道能否优化到$n^2$以下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
// NTT Version
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef long double ldb;
typedef pair <int,int> Pii;
#define reg register
#define pb push_back
#define mp make_pair
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

bool Mbe;

const int N=2050,P=998244353;

int n,m,Z; ll k;
int A[N],B[N];
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

int C[N][N];
int rev[N],T[N],IT[N];
void NTT(int n,int *a,int f){
	rep(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	static int e[N>>1];	
	for(reg int i=e[0]=1;i<n;i<<=1) {
		ll t=f==1?T[i]:IT[i];
		//qpow(f==1?3:(P+1)/3,(P-1)/i/2);
		for(reg int j=i-2;j>=0;j-=2) e[j+1]=t*(e[j]=e[j>>1])%P;
		for(reg int l=0;l<n;l+=i*2) {
			for(reg int j=l;j<l+i;++j) {
				reg int t=1ll*a[j+i]*e[j-l]%P;
				a[j+i]=a[j]-t,Mod2(a[j+i]);
				a[j]+=t,Mod1(a[j]);
			}
		}
	}
	if(f==-1) {
		ll base=qpow(n);
		rep(i,0,n-1) a[i]=a[i]*base%P;
	}
}

int Init(int n){
	int R=1,cc=-1;
	while(R<=n) R<<=1,cc++;
	rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
	return R;
}


int mk1[N],mk2[N];
void FWT(int n,int *A,int *B,int f,int *mk){
		static int X[N],Y[N];
		rep(i,0,n) {
			memset(X,0,sizeof X),memset(Y,0,sizeof Y);
			rep(k,0,i) X[i-k]=(k&1)?P-C[i][k]:C[i][k];
			rep(j,0,n) Y[j]=A[j];
			int R=Init(n+i+1);
			NTT(R,X,1),NTT(R,Y,1);
			rep(j,0,R-1) X[j]=1ll*X[j]*Y[j]%P;
			NTT(R,X,-1);
			rep(j,0,n-i) B[i]=(B[i]+1ll*C[n-i][j]*X[i+j])%P;
		}
}

void Solve(int n,int *A,int *mk) {
	static int X[N],Y[N];
	memset(X,0,sizeof X),memset(Y,0,sizeof Y);
	X[1]=1,FWT(n,X,Y,1,mk);
	rep(i,0,n) Y[i]=qpow(Y[i],k);
	memset(X,0,sizeof X);
	FWT(n,Y,X,2,mk);
	ll base=qpow(qpow(2,n),P-2);
	rep(i,0,n) A[i]=X[i]*base%P*C[n][i]%P;
}


bool Med;
int main(){
	//fprintf(stderr,"%.2lf\n",(&Med-&Mbe)/1024.0/1024.0);
	freopen("clone.in","r",stdin),freopen("clone.out","w",stdout);
	rep(i,0,N-1) rep(j,C[i][0]=1,i) C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
	for(int i=1;i<N;i<<=1) {
		T[i]=qpow(3,(P-1)/i/2);
		IT[i]=qpow((P+1)/3,(P-1)/i/2);
	}
	scanf("%d%d%lld%d",&n,&m,&k,&Z);
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		mk1[i]=1,mk2[j]=1;
	}
	Solve(n,A,mk1),Solve(m,B,mk2);
	int ans=0;
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		ans=(ans+1ll*A[i]*B[j])%P;
	}
	printf("%d\n",ans);
}


// n^2
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

const int N=2050,P=998244353;

int n,m,Z; ll k;
int A[N],B[N];
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int C[N][N],W[N][N],mk1[N],mk2[N];

void Solve(int n,int *A,int *mk) {
	static int X[N],Y[N];
	rep(i,0,n) W[0][i]=C[n][i];
	rep(i,1,n) {
		rep(j,0,n) W[i][j]=W[i-1][j]-(j?W[i][j-1]:0),Mod2(W[i][j]);
		drep(j,n,0) W[i][j]=(j?-W[i][j-1]:0)+W[i][j],Mod2(W[i][j]);
	}

	memset(Y,0,sizeof Y);
	rep(i,0,n) X[i]=qpow(W[i][1],k);
	rep(i,0,n) if(mk[i]) rep(j,0,n) Y[i]=(Y[i]+1ll*W[i][j]*X[j])%P;
	ll base=qpow(qpow(2,n),P-2);
	rep(i,0,n) A[i]=Y[i]*base%P*C[n][i]%P;
}

int main(){
	freopen("clone.in","r",stdin),freopen("clone.out","w",stdout);
	scanf("%d%d%lld%d",&n,&m,&k,&Z);
	rep(i,0,max(n,m)) rep(j,C[i][0]=1,i) C[i][j]=C[i-1][j-1]+C[i-1][j],Mod1(C[i][j]);
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		mk1[i]=mk2[j]=1;
	}
	Solve(n,A,mk1),Solve(m,B,mk2);
	int ans=0;
	rep(i,0,n) rep(j,0,m) {
		if(i*(m-j)+j*(n-i)!=Z) continue;
		ans=(ans+1ll*A[i]*B[j])%P;
	}
	printf("%d\n",ans);
}