Codechef Oct chanllenge Queries on Matrix-JIIT

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 ightarrow dp_{i+1,j-1})

(dp_{i,j}cdot (n-j) ightarrow dp_{i+1,j+1})

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

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

(O(n^3))

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

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

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

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

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

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

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

考虑在这样的多项式上模拟原先的( ext{FWT})过程

---

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

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

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

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

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

(O(n^2log n))

用( ext{NTT})优化上式

(O(n^2))

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

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

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

实际复杂度为(O(n^2+nlog q))

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

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

// 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
",(&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
",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
",ans);
}