【2020五校联考NOIP #2】矩阵

咕咕咕到现在~
题面传送门

题意：
给出一个 \(n\times n\) 的矩阵 \(A\)。要你求有多少个 \(n\times n\) 的矩阵 \(B\) 满足：

• 每一行都是 \(1\) 到 \(n\) 的排列。
• 对于任意 \(1\leq i\lt n\)，\(1\leq j\leq n\)，\(B_{i,j}\neq B_{i+1,j}\)
• 定义 \(f(A)\) 为矩阵 \(A\) 从上到下，从左到右拼接而成的序列，即 \(f(A)_{(i-1)\times n+j}=A_{i,j}\)，那么必须有 \(f(B)\) 的字典序小于 \(f(A)\) 的字典序。
  答案对 \(998244353\) 取模。

我们考虑条件三。“\(f(B)\) 的字典序小于 \(f(A)\) 的字典序”意味着我们可以枚举第一个 \(f(A)_i<f(B)_i\) 的位置然后统计答案。
假设第一个不相等的位置为 \(i\) 行 \(j\) 列。那么矩阵 \(B\) 的前 \(i-1\) 行上的数，以及第 \(i\) 行的前 \(j-1\) 个数，都已经确定下来了，我们的任务就是填好剩下来的 \(n^2-(i-1)\times n-(j-1)\) 个数。
不难注意到，对于后 \(n-i\) 行每行的方案数是固定的，即 \(f_n\)，其中 \(f_i\) 表示长度为 \(i\) 的满足 \(p_j \neq j\) 的排列个数，也就是一个简单的错排数问题。
接下来考虑如何计算填好第 \(i\) 行剩余的 \(n-j+1\) 个数的方案数，我们可以枚举第 \(j\) 个数是什么。那么我们需把这一行剩下来的 \(n-j\) 个数塞入这 \(n-j\) 个位置。再，假设 \(A_{i-1,j+1},A_{i-1,j+2},\dots,A_{i-1,n}\) 分别记作 \(x_1,x_2,\dots,x_{n-j}\)，第 \(i\) 行剩余的数分别记作 \(y_1,y_2,\dots,y_{n-j}\)，要求有多少个 \(1\) 到 \(n-j\) 的排列 \(p\) 满足 \(x_i \neq y_{p_i}\)。
其实方案数与 \(x_i,y_i\) 具体是什么数不重要，我们关心的只是有多少个数在 \(x,y\) 中同时出现，也就是官方题解中的“有限制的位置”与“无限制的位置”。
然后就可以 \(dp\) 了，设 \(dp_{i,j}\) 为有 \(i\) 个有限制的位置和 \(j\) 个无限制的位置的方案数。
若 \(j=0\)，显然 \(dp_{i,j}=f_i\)，直接上公式 \(f_i=\sum\limits_{j=0}^i\binom{i}{j}\times(-1)^{i-j}\)
若 \(j \neq 0\)，那么可以用类似组合数递推的方式。枚举最后一个无限制的位置上填的数，如果填了一个无限制的数，有 \(j\) 种方案，无限制的数变为 \(j-1\)，即 \(j \times dp_{i,j-1}\)；如果填了一个有限制的数，有 \(i\) 种方案，无限制的位置还是 \(j\) 个，原来有限制的数变为无限制的数，即 \(i \times dp_{i-1,j}\)。故 \(dp_{i,j}=j \times dp_{i,j-1}+i \times dp_{i-1,j}\)。
这样做是 \(n^3\) 的。不过可以进一步优化：枚举这个位置填的是有限制的数还是没限制的数，用树状数组/线段树维护选法，乘上一个 \(dp\) 系数即可。

#include <bits/stdc++.h>
using namespace std;
#define fz(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define fill0(a) memset(a,0,sizeof(a))
#define fill1(a) memset(a,-1,sizeof(a))
#define fillbig(a) memset(a,63,sizeof(a))
#define fi first
#define se second
#define int long long
const int MOD=998244353;
int n,a[2005][2005],f[2005],c[2005][2005],fac[2005],dp[2005][2005],b[2005][2005];
bool vis1[2005],vis2[2005];
int qpow(int x,int e){
	int ans=1;
	while(e){
		if(e&1) ans=ans*x%MOD;
		x=x*x%MOD;e>>=1;
	}
	return ans;
}
struct node{
	int l,r,cnt[4],val;
} s[2005<<2];
inline void pushup(int k){
	fz(i,0,3) s[k].cnt[i]=s[k<<1].cnt[i]+s[k<<1|1].cnt[i];
}
inline void build(int k,int l,int r){
//	cout<<k<<" "<<l<<" "<<r<<endl;
	s[k].l=l;s[k].r=r;fill0(s[k].cnt);s[k].val=0;s[k].cnt[0]=r-l+1;
	if(l==r) return;
	int mid=(l+r)>>1;
	build(k<<1,l,mid);build(k<<1|1,mid+1,r);
}
inline void add(int k,int ind,int x){
	if(s[k].l==s[k].r){
		s[k].cnt[s[k].val]--;s[k].val+=x;s[k].cnt[s[k].val]++;
		return;
	}
	int mid=(s[k].l+s[k].r)>>1;
	if(ind<=mid) add(k<<1,ind,x);
	else add(k<<1|1,ind,x);
	pushup(k);
}
inline int query(int k,int l,int r,int x){
	if(l>r||!l||!r) return 0;
	if(l<=s[k].l&&s[k].r<=r) return s[k].cnt[x];
	int mid=(s[k].l+s[k].r)>>1;
	if(r<=mid) return query(k<<1,l,r,x);
	else if(l>mid) return query(k<<1|1,l,r,x);
	else return query(k<<1,l,mid,x)+query(k<<1|1,mid+1,r,x);
}
signed main(){
	scanf("%d",&n);
	fz(i,0,n){
		c[i][0]=1;
		fz(j,1,i) c[i][j]=(c[i-1][j]+c[i-1][j-1])%MOD;
	}
	fac[0]=1;fz(i,1,n) fac[i]=fac[i-1]*i%MOD;
	fz(i,0,n){
		fz(j,0,i){
			if((j&1)==(i&1)) f[i]=(f[i]+c[i][j]*fac[j]%MOD)%MOD;
			else f[i]=(f[i]-c[i][j]*fac[j]%MOD+MOD)%MOD;
		}
//		cout<<i<<" "<<f[i]<<endl;
	}
	for(int i=0;i<=n;i++){
		dp[i][0]=f[i];
		for(int j=1;i+j<=n;j++){
			dp[i][j]=(dp[i][j-1]*j%MOD+dp[i-1][j]*i%MOD)%MOD;
//			printf("%d %d %d\n",i,j,dp[i][j]); 
		}
	}
	fz(i,1,n) fz(j,1,n) scanf("%d",&a[i][j]);
	fz(i,1,n-1){
		fill0(vis1);fill0(vis2);
		int both=0;
		fd(j,n,1){
			vis1[a[i][j]]=1;if(vis2[a[i][j]]) both++;
			vis2[a[i+1][j]]=1;if(vis1[a[i+1][j]]) both++;
			b[i][j]=both;
		}
//		fz(j,1,n) cout<<b[i][j]<<" ";puts("");
	}
	int ans=0;
	fz(i,1,n){
		build(1,1,n);int sum=0;
		fd(j,n,1){
//			cout<<i<<" "<<j<<endl;
			add(1,a[i][j],1);if(i!=1) add(1,a[i-1][j],2);
			int cnt1=query(1,1,a[i][j]-1,1)-((a[i-1][j]<a[i][j])?(query(1,a[i-1][j],a[i-1][j],1)):0);
			int cnt2=query(1,1,a[i][j]-1,3)-((a[i-1][j]<a[i][j])?(query(1,a[i-1][j],a[i-1][j],3)):0);
			int flg=query(1,a[i][j],a[i][j],2)+query(1,a[i][j],a[i][j],3);
			sum=(sum+cnt2*dp[b[i-1][j+1]-1+flg][n-j-(b[i-1][j+1]-1+flg)]%MOD+cnt1*dp[b[i-1][j+1]+flg][n-j-(b[i-1][j+1]+flg)]%MOD)%MOD;
//			cout<<i<<" "<<j<<" "<<cnt1<<" "<<cnt2<<" "<<b[i-1][j]<<endl;
		}
		ans=(ans+sum*qpow(f[n],n-i)%MOD)%MOD;
	}
	printf("%lld\n",ans);
	return 0;
}