Problem:


Score : 1700 points

Problem Statement

You have two strings A = A_1 A_2 ... A_n and B = B_1 B_2 ... B_n of the same length consisting of 0 and 1.
The number of 1's in A and B is equal.

You've decided to transform A using the following algorithm:

Let a_1, a_2, ..., a_k be the indices of 1's in A.

Let b_1, b_2, ..., b_k be the indices of 1's in B.

Replace a and b with their random permutations, chosen independently and uniformly.

For each i from 1 to k, in order, swap A_{a_i} and A_{b_i}.

Let P be the probability that strings A and B become equal after the procedure above.

Let Z = P \times (k!)^2. Clearly, Z is an integer.

Find Z modulo 998244353.

Constraints

1 \leq |A| = |B| \leq 10,000

A and B consist of 0 and 1.

A and B contain the same number of 1's.

A and B contain at least one 1.

Partial Score

1200 points will be awarded for passing the testset satisfying 1 \leq |A| = |B| \leq 500.

Input

Input is given from Standard Input in the following format:

A
B

Output

Print the value of Z modulo 998244353.

Sample Input 1

1010
1100

Sample Output 1

3

After the first two steps, a = [1, 3] and b = [1, 2]. There are 4 possible scenarios after shuffling a and b:

a = [1, 3], b = [1, 2]. Initially, A = 1010. After swap(A_1, A_1), A = 1010. After swap(A_3, A_2), A = 1100.

a = [1, 3], b = [2, 1]. Initially, A = 1010. After swap(A_1, A_2), A = 0110. After swap(A_3, A_1), A = 1100.

a = [3, 1], b = [1, 2]. Initially, A = 1010. After swap(A_3, A_1), A = 1010. After swap(A_1, A_2), A = 0110.

a = [3, 1], b = [2, 1]. Initially, A = 1010. After swap(A_3, A_2), A = 1100. After swap(A_1, A_1), A = 1100.

Out of 4 scenarios, 3 of them result in A = B. Therefore, P = 3 / 4, and Z = 3.

Sample Input 2

01001
01001

Sample Output 2

4

No swap ever changes A, so we'll always have A = B.

Sample Input 3

101010
010101

Sample Output 3

36

Every possible sequence of three swaps puts the 1's in A into the right places.

Sample Input 4

1101011011110
0111101011101

Sample Output 4

932171449



Sample input:
1010
1100


Sample output:
3

Submission status: Accepted

Solution code:
def ilog2(n):
    return 0 if n <= 0 else n.bit_length() - 1
def pack(pack, shamt):
    size = len(pack)
    while size > 1:
        npack = []
        for i in range(0, size - 1, 2):
            npack += [pack[i] | (pack[i+1] << shamt)]
        if size & 1:
            npack += [pack[-1]]
        pack, size, shamt = npack, (size + 1) >> 1, shamt << 1
    return pack[0]
def unpack(M, size, shamt):
    s, sizes = size, []
    while s > 1:
        sizes += [s]
        s = (s + 1) >> 1
    ret = [M]
    for size in sizes[::-1]:
        mask, nret = (1 << shamt) - 1, []
        for c in ret:
            nret += [c & mask, c >> shamt]
        ret, shamt = nret[:size], shamt >> 1
    return ret
def poly_mul_mod(f, g, mod):
    size = min(len(f), len(g))
    shift = ((mod - 1) ** 2 * size).bit_length()
    rsize = len(f) + len(g) - 1
    h = unpack(pack(f, shift) * pack(g, shift), rsize, shift * (1 << ilog2(rsize - 1)))
    return [int(x % mod) for x in h]
mod, N = 998244353, 10 ** 4 + 10
def inv(n): return pow(n, mod - 2, mod)
def poly_power_mod(f, n, mx):
	ret = [1]
	while n:
		if n & 1 : ret = poly_mul_mod(ret, f, mod)[:mx]
		f = poly_mul_mod(f, f, mod)[:mx]
		n >>= 1
	return ret	

fact, ifact, sm, poly = [1] * N, [1] * N, 1, [1] * N
for i in range(1, N): 
	fact[i], ifact[i] = (fact[i - 1] * i) % mod, (ifact[i - 1] * inv(i)) % mod
	poly[i - 1] = ifact[i]

A, B, n, m = raw_input(), raw_input(), 0, 0
for i in range(len(A)):
	if(A[i] == '1'):
		n += 1
		if(B[i] != '1'): m += 1		

ans, poly = 0, poly_power_mod(poly, m, n + 1)
for k in range(n - m + 1):
	term = fact[n - m - k] * fact[n - m - k] % mod
	if k < len(poly):
		term *= poly[k] * fact[m]
	else : term = 0	
	term *= fact[n - m] * fact[n] * ifact[n - m - k] * ifact[n - m - k]
	ans += term

print ans % mod

Submission status: Runtime Error

Solution code:
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<vector>

using namespace std;

#define Whats(x) cout<<#x<<" is "<<x<<endl
#define Divhe() cout<<"<<<<<<<<<<<<<<"<<endl
#define Divshe() cout<<">>>>>>>>>>>>>>>"<<endl
#define Oops() cout<<"!!!!!!!!!!!!!!"<<endl

#define Poly vector<int>

ostream& operator<<(ostream& os,Poly p)
{
	for(int i:p)
		os<<i<<" ";
	return os;
}

const int N=20000;
const int MOD=998244353;
const int g=3;

int Mul(int a,int b)
{
	unsigned long long x=(long long)a*b;
	unsigned xh=(unsigned)(x>>32),xl=(unsigned)x,d,m;
	asm
		(
		 "divl %4;\n\t"
		 :"=a"(d),"=d"(m)
		 :"d"(xh),"a"(xl),"r"(MOD)
		);
	return m;
}

int Qpow(int x,int y=MOD-2)
{
	int res=1;
	for(;y;y>>=1,x=Mul(x,x))if(y&1)
		res=Mul(res,x);
	return res;
}

int U(int x,int y)
{
	return ((x+=y)>=MOD)?(x-MOD):x;
}

void SU(int& x,int y)
{
	((x+=y)>=MOD)?(x-=MOD):0;
}

int D(int x,int y)
{
	return ((x-=y)<0)?(x+MOD):x;
}

void SD(int& x,int y)
{
	((x-=y)<0)?(x+=MOD):0;
}

void FFT(Poly& A,bool fl)
{
	int L=A.size();
	if(fl)
	{
		int t=Qpow(L);
		for(int& i:A)
			i=Mul(i,t);
		reverse(A.begin()+1,A.end());
	}
	for(int i=1,j=L>>1,k;i<L;++i,j^=k)
	{
		if(i<j)
			swap(A[i],A[j]);
		k=L>>1;
		while(j&k)
		{
			j^=k;
			k>>=1;
		}
	}

	static int w[N]={1};

	for(int i=1;i<L;i<<=1)
	{
		int t=Qpow(g,(MOD-1)/(i<<1));
		for(int j=1;j<i;++j)
			w[j]=Mul(w[j-1],t);
		for(int j=0;j<L;j+=(i<<1))
		{
			for(int k=0;k<i;++k)
			{
				t=Mul(A[i+j+k],w[k]);
				A[i+j+k]=D(A[j+k],t);
				SU(A[j+k],t);
			}
		}
	}
}

int n;
char A[N],B[N];

Poly operator*(Poly A,Poly B)
{
	int need=A.size()+B.size()-1,L,q=A.size();
	for(L=1;L<need;L<<=1);
	A.resize(L);
	FFT(A,false);
	B.resize(L);
	FFT(B,false);
	for(int i=0;i<L;++i)
		A[i]=Mul(A[i],B[i]);
	FFT(A,true);
	A.resize(q);
	return A;
}

Poly f;

Poly Qpow(Poly f,int y)
{
	int L=f.size();
	Poly res;
	res.resize(L);
	res[0]=1;
	for(;y;y>>=1,f=f*f)if(y&1)
			res=res*f;
	return res;
}

int jc[N],njc[N];

void Prework()
{
	jc[0]=1;
	for(int i=1;i<=n;++i)
		jc[i]=Mul(jc[i-1],i);
	njc[n]=Qpow(jc[n]);
	for(int i=n-1;i>=0;--i)
		njc[i]=Mul(njc[i+1],i+1);
}

int main()
{
	scanf("%s%s",A,B);
	n=strlen(A);
	int x=0,y=0;
	Prework();
	for(int i=0;i<n;++i)
		if(A[i]=='1')
			(B[i]=='0')?(++x):(++y);
	f.resize(y+1);
	for(int i=0;i<=y;++i)
		f[i]=njc[i+1];
	f=Qpow(f,x);
	int ans=0;
	for(int i=0;i<=y;++i)
		SU(ans,f[i]);
	ans=Mul(ans,Mul(Mul(jc[x],jc[y]),jc[x+y]));
	printf("%d",ans);
	return 0;
}