约定：为方便描述，将 “和为 $2$ 的行”简称为 $2$ 类行，“和为 $1$ 的行”简称为 $1$ 类行，列同理。

设有 $c_1$ 个 $1$ 类行，$c_2$ 个 $2$ 类行，$c_3$ 个 $1$ 类列， $c_4$ 个 $2$ 类列。

基本分析

我们可以固定每行的和为 $R_i$，将 $R_i$ 拆开放在行中的不同位置，来满足每列的和为 $C_i$ 的要求。

因为 $0\leq R_i,C_i \leq 2$，容易想到使用动态规划。

具体实现

状态

设 $f_{i,j,k,l}$ 表示还剩下 $i$ 个 $1$ 类行和 $j$ 个 $2$ 类行，来满足 $k$ 个 $1$ 类列和 $l$ 个 $2$ 类列的方案数。

显然，若 $i+2\times j \ne k+2\times l$，肯定无解，所以所有的 $(i,j,k,l)$ 一定满足 $i+2\times j = k+2\times l$。

转移

对于所有的 $1$ 类行，我们放到最后来考虑。

对于每一个 $2$ 类行，它有这四种转移：

• 拆成 $1+1$ 来满足两个 $1$ 类列，方案数为 $f_{i,j-1,k-2,l}\times \binom{k}{2}$。

• 拆成 $1+1$ 来满足一个 $1$ 类列，并将一个 $2$ 类列变成 $1$ 类列，方案数为 $f_{i,j-1,k,l-1}\times k\times l$。

• 拆成 $1+1$ 将两个 $2$ 类列变成 $1$ 类列，方案数为 $f_{i,j-1,k,l-2}\times \binom{l}{2}$。

• 满足一个 $2$ 类列，方案数为 $f_{i,j-1,k,l-1}\times l$。

边界条件

$j=0$ 时，剩下 $i$ 个 $1$ 类行，满足 $k$ 个 $1$ 类列和 $l$ 个 $2$ 类列。

即 $f_{i,0,k,l}=\binom{k}{1}\binom{k-1}{1}\dots\binom{2\times l+1}{1}\binom{2\times l}{2}\binom{2\times l-2}{2}\dots\binom{2}{2}=\frac{k!}{2^l}$。

时间复杂度

状态数为 $O(n^4)$，转移 $O(1)$，时间复杂度为 $O(n^4)$，无法通过。

优化状态数

所有状态的第一维 $i$ 均是不变的，可以省掉这一维。

于是现在的状态 $f_{j,k,l}$ 表示还剩下 $c_1$ 个 $1$ 类行和 $j$ 个 $2$ 类行，来满足 $k$ 个 $1$ 类列和 $l$ 个 $2$ 类列的方案数。

还记得在这篇题解最开始提到的那个等式吗： $i+2\times j = k+2\times l$。

所以 $k$ 可以表示为 $i+2\times j -2\times l$。

那么现在的状态即 $f_{j,l}$ 表示还剩下 $c_1$ 个 $1$ 类行和 $j$ 个 $2$ 类行，来满足 $c_1+2\times j -2\times l$ 个 $1$ 类列和 $l$ 个 $2$ 类列的方案数。

转移

设 $k=i+2\times j -2\times l$。

于是 $f_{j,l}=f_{j-1,l}\times \binom{k}{2} + f_{j-1,l-1}\times k\times l+ f_{j-1,l-2}\times \binom{l}{2} + f_{j-1,l-1}\times l$。

边界条件

$f_{0,l}=\frac{c_1!}{2^l}$。

时间复杂度

状态数为 $O(n^2)$，转移 $O(1)$，时间复杂度为 $O(n^2)$，可以通过。

代码（请勿抄袭）

#include<bits/stdc++.h>
#define mod 998244353
#define int long long
using namespace std;
namespace Binom//组合数板子
{
	#define N 50005//注意预处理需要处理2倍n长度
	int st[N],qp[N];
	int qpow(int x,int n)
	{
		if(x<0) return 0;
		if(!n) return 1;
		int t=qpow(x,n/2);
		if(n&1) return t*t%mod*x%mod;
		return t*t%mod;
	}
	void C_init()
	{
		st[0]=1;
		for(int i=1;i<N;i++) st[i]=st[i-1]*i%mod;
		qp[N-1]=qpow(st[N-1],mod-2);
		for(int i=N-2;i>=0;i--) qp[i]=qp[i+1]*(i+1)%mod;
	}
	int C(int n,int m)
	{
		if(n<0||m<0||m>n) return 0;
		return st[n]*qp[n-m]%mod*qp[m]%mod;
	}
	int A(int n,int m)
	{
		if(n<0||m<0||m>n) return 0;
		return st[n]*qp[n-m]%mod;
	}
	#undef N
}
using namespace Binom;
#define N 5005
int n,r[N],c[N],f[N][N];
void add(int &x,int y){x=(x+y)%mod;}
signed main()
{
	C_init();
	scanf("%lld",&n);
	for(int i=1;i<=n;i++) scanf("%lld",&r[i]);
	for(int i=1;i<=n;i++) scanf("%lld",&c[i]);
	int c1=0,c2=0,c3=0,c4=0;
	for(int i=1;i<=n;i++)
	{
		if(r[i]==1) c1++;
		if(r[i]==2) c2++;
		if(c[i]==1) c3++;
		if(c[i]==2) c4++;
	}
	if(c1+c2*2!=c3+c4*2) return puts("0"),0;
	for(int i=0;i<=c4;i++) f[0][i]=st[c1]*qpow(qpow(2,i),mod-2)%mod;//边界条件
	for(int j=1;j<=c2;j++)
		for(int l=0;l<=c4;l++)//转移
		{
			int k=c1+2*j-2*l;
			if(k>=2) add(f[j][l],f[j-1][l]*C(k,2)%mod);//拆1+1给1+1
			if(l&&k) add(f[j][l],f[j-1][l-1]*k%mod*l%mod);//拆1+1给1+2
			if(l>=2) add(f[j][l],f[j-1][l-2]*C(l,2)%mod);//拆1+1给2+2
			if(l) add(f[j][l],f[j-1][l-1]*l%mod);//给2
		}
	printf("%lld",f[c2][c4]);
	return 0;
}

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=998244353;
ll n,a[5010],b[5010],f[5010][5010],ff[10010];
map<pair<ll,ll>,ll> mp,mp2;
ll qpow(ll a,ll b){
	ll ans=1;
	for(;b;b>>=1,a=a*a%mod)if(b&1)ans=ans*a%mod;
	return ans;
}
int main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cin>>n;
	ff[0]=1;
	for(int i=1;i<=1e4;i++){
		ff[i]=ff[i-1]*i%mod;
	}
	ll a2=0,a1=0,b2=0,b1=0,as=0,bs=0;
	for(int i=1;i<=n;i++){
		cin>>a[i];as+=a[i];
		if(a[i]==1)a1++;
		if(a[i]==2)a2++;
	}
	for(int i=1;i<=n;i++){
		cin>>b[i];bs+=b[i];
		if(b[i]==1)b1++;
		if(b[i]==2)b2++;
	}
	if(as!=bs){
		cout<<0;
		return 0;
	}
	if(a2&&!b2){
		swap(a1,b1);
		swap(a2,b2);
	}
	ll ans=0;
	mp[{b1,b2}]=1;
	if(a2){
		while(a2--){
			mp2.clear();
			for(auto i:mp){
				ll x=i.first.first,y=i.first.second,v=i.second;
				if(x>=2)mp2[{x-2,y}]=(mp2[{x-2,y}]+v*((x*(x-1)/2)%mod)%mod)%mod;
				if(y>0)mp2[{x,y-1}]=(mp2[{x,y-1}]+v*(x+1)%mod*y%mod)%mod;
				if(y>1)mp2[{x+2,y-2}]=(mp2[{x+2,y-2}]+v*((y*(y-1)/2)%mod)%mod)%mod;
			}
			mp.clear();
			for(auto i:mp2){
				ll x=i.first.first,y=i.first.second,v=i.second;
				mp[{x,y}]=v;
			}
		}
	}
	for(auto i:mp){
		ll x=i.first.first,y=i.first.second,v=i.second;
		ans=(ans+v*ff[x+y*2]%mod*qpow(qpow(2,y),mod-2)%mod)%mod;
	}
	cout<<ans;
	return 0;
}