代码

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
#define rep(i, j, k) for(int i = (j); i <= (k); i++)
#define per(i, j, k) for(int i = (j); i >= (k); i--)
#define pb emplace_back
#define fi first
#define se second
using vi = vector<int>;
using pi = pair<int, int>;

template<typename T0, typename T1> bool chmin(T0 &x, const T1 &y){
	if(y < x){x = y; return true;} return false;
}
template<typename T0, typename T1> bool chmax(T0 &x, const T1 &y){
	if(x < y){x = y; return true;} return false;
}

template<typename T> void debug(char *s, T x){
	cerr << s <<" = "<< x <<endl;
}
template<typename T, typename ...Ar> void debug(char *s, T x, Ar... y){
	int dep = 0;
	while(!(*s == ',' && dep == 0)){
		if(*s == '(') dep++;
		if(*s == ')') dep--;
		cerr << *s++;
	}
	cerr <<" = "<< x <<",";
	debug(s + 1, y...);
}
#define gdb(...) debug((char*)#__VA_ARGS__, __VA_ARGS__)

using u32 = uint32_t;
using u64 = uint64_t;
constexpr int mod = 998244353;
struct mint{
	u32 x;
	
	mint(): x(0){}
	mint(int _x){
		_x %= mod;
		if(_x < 0) _x += mod;
		x = _x;
	}

	u32 val()const {
		return x;
	}
	mint qpow(int y = mod - 2)const {
		assert(y >= 0);
		mint t = *this, res = 1;
		while(y){
			if(y%2) res *= t;
			t *= t;
			y /= 2;
		}
		return res;
	}

	mint& operator += (const mint &B){
		if((x += B.x) >= mod) x -= mod;
		return *this;
	}
	mint& operator -= (const mint &B){
		if((x -= B.x) >= mod) x += mod;
		return *this;
	}
	mint& operator *= (const mint &B){
		x = (u64)x * B.x % mod;
		return *this;
	}
	mint& operator /= (const mint &B){
		return *this *= B.qpow();
	}
	friend mint operator + (const mint &A, const mint &B){
		return mint(A) += B;
	}
	friend mint operator - (const mint &A, const mint &B){
		return mint(A) -= B;
	}
	friend mint operator * (const mint &A, const mint &B){
		return mint(A) *= B;
	}
	friend mint operator / (const mint &A, const mint &B){
		return mint(A) /= B;
	}
	mint operator - ()const {
		return mint() - *this;
	}
};

constexpr int LG = 21, S = 1 << LG;
vector<mint> P, fac, ifac, iv, ip2;

auto INIT = [](){
	P.resize(S);
	rep(i, 0, LG - 1){
		mint stp = mint(3).qpow(mod >> (i + 1));
		P[1 << i] = 1;
		rep(j, (1 << i) + 1, (2 << i) - 1){
			P[j] = P[j - 1] * stp;
		}
	}

	fac.resize(S);
	ifac.resize(S);
	iv.resize(S);
	fac[0] = iv[0] = 1;
	rep(i, 1, S - 1){
		fac[i] = fac[i - 1] * i;
	}
	ifac[S - 1] = 1 / fac[S - 1];
	per(i, S - 1, 1){
		ifac[i - 1] = ifac[i] * i;
		iv[i] = ifac[i] * fac[i - 1];
	}

	ip2.resize(LG + 1);
	ip2[0] = 1;
	rep(i, 1, LG){
		ip2[i] = ip2[i - 1] / 2;
	}
	
	return true;
}();

struct poly : vector<mint>{
	using vector<mint>::vector;
	
	void DIF(int lg){
		for(int i = 1 << (lg - 1); i ; i /= 2){
			for(int j = 0; j < (1 << lg); j += 2 * i){
				rep(k, 0, i - 1){
					mint A = (*this)[j + k], B = (*this)[j + k + i];
					(*this)[j + k] = A + B, (*this)[j + k + i] = (A - B) * P[i + k];
				}
			}
		}
	}

	void DIT(int lg){
		for(int i = 1; i < (1 << lg); i *= 2){
			for(int j = 0; j < (1 << lg); j += 2 * i){
				rep(k, 0, i - 1){
					mint A = (*this)[j + k], B = (*this)[j + k + i] * P[i + k];
					(*this)[j + k] = A + B, (*this)[j + k + i] = A - B;
				}
			}
		}
		reverse(this -> begin() + 1, this -> end());
	}

	friend poly der(const poly &A){
		poly res = A;
		rep(i, 1, (int)res.size() - 1){
			res[i - 1] = res[i] * i;
		}
		res.pop_back();
		return res;
	}

	friend poly ints(const poly &A){
		poly res = A;
		res.pb();
		per(i, (int)res.size() - 1, 0){
			res[i] = res[i - 1] * iv[i];
		}
		res[0] = 0;
		return res;
	}

	poly& operator += (const poly &B){
		if(B.size() > this -> size()){
			this -> resize(B.size());
		}
		rep(i, 0, (int)B.size() - 1){
			(*this)[i] += B[i];
		}
		return *this;
	}

	friend poly operator + (const poly &A, const poly &B){
		return poly(A) += B;
	}

	poly& operator -= (const poly &B){
		if(B.size() > this -> size()){
			this -> resize(B.size());
		}
		rep(i, 0, (int)B.size() - 1){
			(*this)[i] -= B[i];
		}
		return *this;
	}

	friend poly operator - (const poly &A, const poly &B){
		return poly(A) -= B;
	}

	friend poly operator * (const mint &k, const poly &A){
		poly res = A;
		for(auto &x:res) x *= k;
		return res;
	}

	poly& operator *= (const poly &B){
		auto sz = this -> size();
		*this = conv(*this, B);
		this -> resize(sz);
		return *this;
	}

	friend poly operator * (const poly &A, const poly &B){
		return poly(A) *= B;
	}

	poly& operator /= (const poly &B){
		auto sz = this -> size();
		*this = conv(*this, inv(B));
		this -> resize(sz);
		return *this;
	}

	friend poly operator / (const poly &A, const poly &B){
		return poly(A) /= B;
	}

	friend poly conv(const poly &_A, const poly &_B){
		if(_A.size() <= 64 || _B.size() <= 64){
			poly res(_A.size() + _B.size() - 1);
			rep(i, 0, (int)_A.size() - 1){
				rep(j, 0, (int)_B.size() - 1){
					res[i + j] += _A[i] * _B[j];
				}
			}
			return res;
		}

		poly A = _A, B = _B;
		int len = A.size() + B.size() - 1;
		int lg = __lg(len * 2 - 1);
		A.resize(1 << lg), B.resize(1 << lg);
		A.DIF(lg), B.DIF(lg);
		rep(i, 0, (1 << lg) - 1){
			A[i] *= B[i];
		}
		A.DIT(lg);
		A.resize(len);
		for(auto &x:A) x *= ip2[lg];
		return A;
	}

	friend poly inv(const poly &A){
		assert(A[0].val() != 0);
		poly res{mint(1) / A[0]};
		while(res.size() < A.size()){
			int n = res.size(), lg = __lg(n * 4);
			poly buf(1 << lg);
			copy(res.begin(), res.end(), buf.begin());
			buf.DIF(lg);

			poly tmp{A.begin(), A.begin() + min(n * 2, (int)A.size())};
			tmp.resize(n * 4);
			tmp.DIF(lg);
			rep(i, 0, n * 4 - 1) tmp[i] *= buf[i] * buf[i];
			tmp.DIT(lg);

			res.resize(n * 2);
			rep(i, 0, n * 2 - 1){
				res[i] = 2 * res[i] - tmp[i] * ip2[lg];
			}
		}
		res.resize(A.size());
		return res;
	}


	// O(n \log^2 n)
	friend poly exp(const poly &A){
		assert(A[0].val() == 0);
		int n = A.size();
		int lg = __lg(n * 2 - 1);

		poly a = A, res(A.size());
		rep(i, 0, (int)a.size() - 1) a[i] *= i;
		a.resize(1 << lg);

		vector<poly> pre(lg + 1);
		rep(i, 8, lg){
			pre[i].insert(pre[i].end(), a.begin(), a.begin() + (1 << i));
			pre[i].DIF(i);
		}

		res[0] = 1;
		auto slv = [&](auto &self, int l, int r)->void {
			if(r - l <= 128){
				rep(i, l, r - 1){
					rep(j, l, i - 1){
						res[i] += res[j] * a[i - j];
					}
					if(i > 0) res[i] *= iv[i];
				}
				return;
			}
			int mid = (l + r) / 2;
			self(self, l, mid);

			int _lg = __lg((r - l) * 2 - 1);
			poly tmp(res.begin() + l, res.begin() + mid);
			tmp.resize(1 << _lg);
			tmp.DIF(_lg);
			rep(i, 0, (1 << _lg) - 1) tmp[i] *= pre[_lg][i];
			tmp.DIT(_lg);
			rep(i, mid, r - 1) res[i] += tmp[i - l] * ip2[_lg];
			self(self, mid, r);
		};
		slv(slv, 0, n);
		return res;
	}
	
	friend poly ln(const poly &A){
		assert(A[0].val() == 1);
		poly res = der(A) / A;
		return ints(res);
	}
	
	friend poly sqrt(const poly &A){
		assert(A[0].val() == 1);
		return exp(ip2[1] * ln(A));
	}

	friend poly qpow(const poly &A, mint k){
		if(k.val() == 0){
			poly res(A.size());
			res[0] = 1;
			return res;
		}

		auto it = A.begin();
		while(it != A.end() && it -> val() == 0){
			it ++;
		}
		int len = (it - A.begin()) * k.val();
		if(len >= (int)A.size()){
			return poly(A.size());
		} else{
			poly a{it, it + (A.size() - len)};
			mint t = a[0].qpow(k.val());
			a = t * exp(k * ln(1 / a[0] * a));
			a.insert(a.begin(), len, mint());
			return a;
		}
	}
};

namespace SPS{
	using sps = vector<mint>;
	constexpr int N = 21;
	mint fac[N], ifac[N], inv[N];
	void init(){
		fac[0] = 1;
		rep(i, 1, N - 1) fac[i] = fac[i - 1] * i;
		ifac[N - 1] = 1 / fac[N - 1];
		per(i, N - 1, 1) ifac[i - 1] = ifac[i] * i, inv[i] = ifac[i] * fac[i - 1];
	}

	void fmt(sps &a, int n){
		for(int i = 1; i < (1<<n); i *= 2){
			for(int j = 0; j < (1<<n); j += 2 * i){
				rep(k, 0, i - 1){
					a[i + j + k] += a[j + k];
				}
			}
		}
	}
	void ifmt(sps &a, int n){
		for(int i = 1; i < (1<<n); i *= 2){
			for(int j = 0; j < (1<<n); j += 2 * i){
				rep(k, 0, i - 1){
					a[i + j + k] -= a[j + k];
				}
			}
		}
	}
	sps sps_conv(const sps &a, const sps &b){
		int n = __lg(a.size());
		vector<sps> A(n + 1, sps(1 << n)), B = A, C = A;
		rep(s, 0, (1 << n) - 1){
			int pc = __builtin_popcount(s);
			A[pc][s] = a[s];
			B[pc][s] = b[s];
		}
		rep(i, 0, n){
			fmt(A[i], n);
			fmt(B[i], n);
		}
		rep(i, 0, n){
			rep(j, 0, i){
				rep(k, 0, (1 << n) - 1){
					C[i][k] += A[j][k] * B[i - j][k];
				}
			}
		}
		rep(i, 0, n){
			ifmt(C[i], n);
		}
		sps res(1 << n);
		rep(i, 0, (1 << n) - 1){
			res[i] = C[__builtin_popcount(i)][i];
		}
		return res;
	}
	sps sps_exp(const sps &a){
		int n = __lg(a.size());
		vector< sps > b(n + 1, sps(1 << n)), c = b;
		rep(i, 0, (1 << n) - 1){
			b[__builtin_popcount(i)][i] = a[i];
		}
		rep(i, 0, n){
			fmt(b[i], n);
		}
		rep(i, 0, n){
			rep(j, 0, (1 << n) - 1){
				b[i][j] *= i;
			}
		}
		fill(c[0].begin(), c[0].end(), 1);
		rep(i, 1, n){
			rep(j, 1, i){
				rep(k, 0, (1 << n) - 1){
					c[i][k] += b[j][k] * c[i - j][k];
				}
			}
			rep(k, 0, (1 << n) - 1){
				c[i][k] *= inv[i];
			}
		}
		rep(i, 0, n){
			ifmt(c[i], n);
		}
		sps res(1 << n);
		rep(i, 0, (1 << n) - 1){
			res[i] = c[__builtin_popcount(i)][i];
		}
		return res;
	}
	sps sps_ln(const sps &a){
		int n = __lg(a.size());
		vector< sps > b(n + 1, sps(1 << n)), c = b;
		rep(i, 0, (1 << n) - 1){
			b[__builtin_popcount(i)][i] = a[i];
		}
		rep(i, 0, n){
			fmt(b[i], n);
		}
		rep(i, 1, n){
			rep(j, 1, i - 1){
				rep(k, 0, (1 << n) - 1){
					c[i][k] += j * c[j][k] * b[i - j][k];
				}
			}
			rep(k, 0, (1 << n) - 1){
				c[i][k] = b[i][k] - c[i][k] * inv[i];
			}
		}
		rep(i, 0, n){
			ifmt(c[i], n);
		}
		sps res(1 << n);
		rep(i, 0, (1 << n) - 1){
			res[i] = c[__builtin_popcount(i)][i];
		}
		return res;
	}
	sps sps_inv(const sps &a){
		int n = __lg(a.size());
		vector< sps > b(n + 1, sps(1 << n)), c = b;
		rep(i, 0, (1 << n) - 1){
			b[__builtin_popcount(i)][i] = a[i];
		}
		rep(i, 0, n){
			fmt(b[i], n);
		}		
		rep(k, 0, (1 << n) - 1){
			c[0][k] = 1 / b[0][k];
		}
		rep(i, 1, n){
			rep(j, 0, i - 1){
				rep(k, 0, (1 << n) - 1){
					c[i][k] -= c[j][k] * b[i - j][k];
				}
			}
			rep(k, 0, (1 << n) - 1){
				c[i][k] *= c[0][k];
			}
		}
		rep(i, 0, n){
			ifmt(c[i], n);
		}
		sps res(1 << n);
		rep(i, 0, (1 << n) - 1){
			res[i] = c[__builtin_popcount(i)][i];
		}
		return res;
	}
	sps sps_comp(const sps &a, const poly &b){
		int n = __lg(a.size());

		vector< vector<mint> > f(n + 1, vector<mint>(1, 1));
		mint fac = 1;
		rep(i, 0, n){
			f[i][0] = b[i] * fac;
			fac *= (i + 1);
		}

		rep(i, 1, n){
			vector< vector<mint> > tmp(n + 1 - i, vector<mint>(1 << i));
			sps F(a.begin() + (1 << (i - 1)), a.begin() + (1 << i));
			rep(j, 0, n - i){
				sps G = sps_conv(F, f[j + 1]);
				copy(f[j].begin(), f[j].end(), tmp[j].begin());
				rep(k, 0, (1 << (i - 1)) - 1){
					tmp[j][k + (1 << (i - 1))] = G[k];
				}
			}
			f = move(tmp);
		}
		return f[0];
	}
}
using sps = vector<mint>;
using SPS::sps_conv;
using SPS::sps_exp;
using SPS::sps_ln;
using SPS::sps_inv;
using SPS::sps_comp;

mint C(int m, int n){
	return fac[m] * ifac[n] * ifac[m - n];
}

void sort(unsigned *a, int n){
	for(int i = 1; i < n; i += 2){
		a[i] = min(a[i] + a[i - 1], a[i] + a[i - 1] - mod);
	}
}

评测结果