Language: C++

#define _USE_MATH_DEFINES
#include <cstdio>
#include <cstring>
#include <cmath>
#include <complex>
using namespace std;

typedef long long ll;

const int MAXN = 10005;
const int MOD = 998244353;

namespace FFT {
  const int MAX = 1 << 20;

	#define REP(i, n) for (int i = 0; i < n; i++)
  typedef long long llint;
  typedef llint value;
  typedef complex<double> comp;

  int N;
  comp omega[MAX];
  comp a1[MAX], a2[MAX];
  comp z1[MAX], z2[MAX];

  void fft(comp *a, comp *z, int m = N) {
    if (m == 1) {
      z[0] = a[0];
    } else {
      int s = N/m;
      m /= 2;

      fft(a, z, m);
      fft(a+s, z+m, m);

      REP(i, m) {
        comp c = omega[s*i] * z[m+i];
        z[m+i] = z[i] - c;
        z[i] += c;
      }
    }
  }
  
  void mult(value *a, value *b, value *c, int len) {
    N = 2*len;
    while (N & (N-1)) ++N;
    assert(N <= MAX);

    REP(i, N) a1[i] = 0;
    REP(i, N) a2[i] = 0;
    REP(i, len) a1[i] = a[i];
    REP(i, len) a2[i] = b[i];

    REP(i, N) omega[i] = polar(1.0, 2*M_PI/N*i);
    fft(a1, z1, N);
    fft(a2, z2, N);

    REP(i, N) omega[i] = comp(1, 0) / omega[i];
    REP(i, N) a1[i] = z1[i] * z2[i] / comp(N, 0);
    fft(a1, z1, N);

    REP(i, 2*len) c[i] = round(z1[i].real());
  }

  void mult_mod(llint *a, llint *b, llint *c, int len, int mod) {
    static llint a0[MAX], a1[MAX];
    static llint b0[MAX], b1[MAX];
    static llint c0[MAX], c1[MAX], c2[MAX];

    REP(i, len) a0[i] = a[i] & 0xFFFF;
    REP(i, len) a1[i] = a[i] >> 16;

    REP(i, len) b0[i] = b[i] & 0xFFFF;
    REP(i, len) b1[i] = b[i] >> 16;

    FFT::mult(a0, b0, c0, len);
    FFT::mult(a1, b1, c2, len);

    REP(i, len) a0[i] += a1[i];
    REP(i, len) b0[i] += b1[i];
    FFT::mult(a0, b0, c1, len);
    REP(i, 2*len) c1[i] -= c0[i] + c2[i];

    REP(i, 2*len) c1[i] %= mod;
    REP(i, 2*len) c2[i] %= mod;
    REP(i, 2*len) c[i] = (c0[i] + (c1[i] << 16) + (c2[i] << 32)) % mod;
  }
}

int N;
char a[MAXN], b[MAXN];
int fact[MAXN], inv[MAXN];
ll base[MAXN], res[MAXN];
 
inline int add(int x, int y) {
	x += y;
	if (x >= MOD)
		x -= MOD;
	return x;
}

inline int mul(int x, int y) {
	return (ll)x * y % MOD;
}

int pot(int x, int y) {
	int res = 1;
	for (; y; y /= 2) {
		if (y % 2)
			res = mul(res, x);
		x = mul(x, x);
	}
	return res;
}

void load() {
	scanf("%s%s", a, b);
}

int solve() {
	N = strlen(a);
	
	fact[0] = inv[0] = 1;
	for (int i = 1; i <= N; i++) {
		fact[i] = mul(fact[i - 1], i);
		inv[i] = pot(fact[i], MOD - 2);
	}
	
	int diff = 0, both = 0;
	for (int i = 0; i < N; i++) {
		diff += a[i] == '1' && b[i] == '0';
		both += a[i] == '1' && b[i] == '1';
	}
	
	for (int i = 0; i <= both; i++)
		base[i] = inv[i + 1];
	res[0] = 1;
	
	for (int i = diff; i; i /= 2) {
		if (i % 2)
			FFT::mult_mod(res, base, res, both + 1, MOD);
		FFT::mult_mod(base, base, base, both + 1, MOD);
	}
	
	int sol = 0;
	for (int i = 0; i <= both; i++)
		sol = add(sol, res[i]);
	
	return mul(sol, mul(mul(fact[both], fact[diff]), fact[diff + both]));
}

int main() {
	load();
	printf("%d\n", solve());
	return 0;
}

Compiler output:
a.cc: In function 'void FFT::mult(value*, value*, value*, int)':
a.cc:47:5: error: 'assert' was not declared in this scope
   47 |     assert(N <= MAX);
      |     ^~~~~~
a.cc:6:1: note: 'assert' is defined in header '<cassert>'; this is probably fixable by adding '#include <cassert>'
    5 | #include <complex>
  +++ |+#include <cassert>
    6 | using namespace std;