Let us define \(f(x)\) by

\[ f(x)= \Bigl(1+A_1x+A_1^2x^2+\cdots+A_1^{B_1}x^{B_1}\Bigr) \Bigl(1+A_2x+A_2^2x^2+\cdots+A_2^{B_2}x^{B_2}\Bigr) \cdots \Bigl(1+A_Nx+A_N^2x^2+\cdots+A_N^{B_N}x^{B_N}\Bigr). \]

Then the desired value is the coefficient of \(x^m\) in \(f(x)\), modulo \(998244353\). Since it is difficult to derive the coefficient directly, let us transform the equation as

\[ f(x)=\frac{\bigl(1-(A_1x)^{B_1+1}\bigr)\bigl(1-(A_2x)^{B_2+1}\bigr)\cdots\bigl(1-(A_Nx)^{B_N+1}\bigr)}{(1-A_1x)(1-A_2x)\cdots (1-A_Nx)} \]

and consider finding the coefficients \(c_1\), \(\ldots\), \(c_N\) satisfying

\[ \frac{1}{(1-A_1x)(1-A_2x)\cdots (1-A_Nx)} =\frac{c_1}{1-A_1x}+ \frac{c_2}{1-A_2x}+\cdots+ \frac{c_N}{1-A_Nx}. \]

By reducing the denominator, we have

\[ 1 = c_1(1-A_2x)(1-A_3x)\cdots (1-A_Nx) +c_2(1-A_1x)(1-A_3x)\cdots (1-A_Nx) +\cdots +c_N(1-A_1x)(1-A_2x)\cdots (1-A_{N-1}x). \]

Here, consider assigning \(x=A_1^{-1},A_2^{-1},\ldots,A_N^{-1}\) to both hand sides. When \(x=A_i^{-1}\), all the terms in the right hand side containing \((1-A_ix)\) evaluates to \(0\), so \(c_i=\Bigl ((1-A_1A_i^{-1})\cdots(1-A_{i-1}A_i^{-1}) (1-A_{i+1}A_i^{-1})\cdots (1-A_NA_i^{-1}) \Bigr)^{-1}\) makes the both hand sides equal at \(x={A_i}^{-1}\). Here, note that \(c_i\) is a finite rational number, since \(A_i\) are pairwise distinct. Now, for such coefficients \(c_1,\ldots, c_N\), \((\text{left hand side})-(\text{right hand side})\) is a polynomial of degree at most \((N-1)\) which evaluates to \(0\) at \(N\) distinct points \(x=A_1^{-1},A_2^{-1},\ldots,A_N^{-1}\), so both hand sides correspond as polynomials. Here, since each \(A_i\) is between \(1\) (inclusive) and \(998244353\) (exclusive), it has an inverse \(A_i^{-1}\) in \(\bmod{ 998244353}\), and the numerator and denominator of \(c_i\) cannot be \(0\) in \(\bmod{998244353}\), so \(c_i\) can be found in \(\bmod{998244353}\). Then,

\[ 1 \equiv c_1(1-A_2x)(1-A_3x)\cdots (1-A_Nx) +c_2(1-A_1x)(1-A_3x)\cdots (1-A_Nx) +\cdots +c_N(1-A_1x)(1-A_2x)\cdots (1-A_{N-1}x) \pmod{998244353}. \]

Here, we can write as \( f(x) \equiv \bigl(1-(A_1x)^{B_1+1}\bigr)\bigl(1-(A_2x)^{B_2+1}\bigr)\cdots\bigl(1-(A_Nx)^{B_N+1}\bigr) \bigl( c_1(1+A_1x+(A_1x)^2+\cdots) +c_2(1+A_2x+(A_2x)^2+\cdots) +\cdots +c_N(1+A_Nx+(A_Nx)^2+\cdots) \bigr) \pmod{998244353} \)

The first half

\[ \bigl(1-(A_1x)^{B_1+1}\bigr)\bigl(1-(A_2x)^{B_2+1}\bigr)\cdots\bigl(1-(A_Nx)^{B_N+1}\bigr) \]

consists of at most \(2^N\) terms. For one of these terms \(d_kx^k\), if \(k > M\) then it cannot contribute to the desired value as a result of multiplication with the second half. Otherwise, the coefficient of \(x^M\) in \(d_kx^k\bigl( c_1(1+A_1x+(A_1x)^2+\cdots) +c_2(1+A_2x+(A_2x)^2+\cdots) +\cdots +c_N(1+A_Nx+(A_Nx)^2+\cdots) \bigr)\) can be found as \(d_k\bigl( c_1A_1^{M-k}+c_2A_2^{M-k}\cdots +c_NA_N^{M-k} \bigr)\). By computing this for each of \(2^N\) terms, the time complexity is \(O(2^N\cdot N\log M)\); for the other part, computing \(d_i\) requires \(O(2^N)\) time and computing \(c_i\) requires \(O(N\log P)\) (where \(P=998244353\) this time), so it will fit in the time limit. Therefore, the problem has been solved. Be careful of the overflows in the implementation.

If you don’t want to treat infinite terms, you can express \(f(x)\) with \(c_i\) as \(\Sigma_{i=1}^n c_i\bigl(1+A_ix+(A_ix)^2+\cdots +(A_ix)^{B_i}) \bigl(1-(A_1x)^{B_1+1}\bigr)\bigl(1-(A_2x)^{B_2+1}\bigr)\cdots \bigl(1-(A_{i-1}x)^{B_{i-1}+1}\bigr)\bigl(1-(A_{i+1}x)^{B_{i+1}+1}\bigr) \cdots\bigl(1-(A_Nx)^{B_N+1}\bigr)\) so that the answer can be found similarly in a total of \(O(N\cdot 2^N\log M)\) time.

Sample code in c++:

#include <bits/stdc++.h>
using namespace std;

#define N 16
#define MOD 998244353
#define ll long long
#define rep(i, n) for(int i = 0; i < n; ++i)

ll modpow(ll x, ll k) {
	ll re = 1;
	while (k > 0) {
		if (k & 1)re = (re*x) % MOD;
		k = (k >> 1);
		x = (x*x) % MOD;
	}
	return re;
}

ll rev(ll x) { return modpow(x, MOD - 2); }

int main() {
	int n;
	ll m;
	ll a[N], b[N];
	ll c[N];
	ll d[(1 << N)], p[(1 << N)];
	ll x, y;
	ll ans = 0;

	cin >> n >> m;
	d[0] = 1LL;
	p[0] = 0LL;

	rep(i, n) {
		cin >> a[i] >> b[i];
		b[i]++;
		x = (MOD - modpow(a[i], b[i])) % MOD;
		rep(j, (1 << i)) {
			d[(1 << i) + j] = (d[j] * x) % MOD;
			p[(1 << i) + j] = p[j] + b[i];
		}
	}

	rep(i, n) {
		x = rev(a[i]);
		c[i] = 1LL;
		rep(j, n) {
			if (j != i) {
				y = (a[j] * x) % MOD;
				c[i] = (c[i] * (MOD + 1 - y)) % MOD;
			}
		}
		c[i] = rev(c[i]);
	}

	rep(i, (1 << n)) {
		if (p[i] <= m) {
			x = 0;
			rep(j, n) {
				x = (x + (modpow(a[j], m - p[i])*c[j])) % MOD;
			}
			ans = (ans + (x*d[i])) % MOD;
		}
	}

	cout << ans << endl;
	return 0;
}