[容斥][组合数学][分治FFT]青春猪头少年不会遇到兔女郎学姐

• source：集训队作业 2019 by 美傲华 mayaohua

一、如果是链而不是环，并且每种颜色分成的段数已经确定

• 考虑如果是链而不是环，并且已经确定第 \(i\) 种颜色在链上恰好组成 \(a_i\) 段，求所有 \(\sum_{i=1}^na_i\) 段排列的方案数

• 也就是对所有的段进行排列（颜色相同的段之间没有顺序），使得任意相邻两段不同色的方案数

• 考虑容斥，即强制所有 \(\sum_{i=1}^n(a_i-1)\) 个分裂点中的部分点被连接起来

• 设 \(f_{i,j}=(-1)^{a_i-j}\binom{a_i-1}{j-1}\)，也就是 \(a_i-1\) 个分裂点有 \(a_i-j\) 个被强制连接，连成 \(j\) 段的方案数

• 则令 EGF \(F_i(x)=\sum_{j}\frac{f_{i,j}}{j!}x^j\)，那么方案数就是 \(\sum_ii![x^i]\prod_{j=1}^nF_j(x)\)

• 分治 NTT 即可，\(O(m\log m\log n)\)，\(m=\sum_{i=1}^na_i\)

二、每种颜色分成的段任意

• 还是链，但是每种颜色分成的段数和段长是任意的，求所有方案的段长乘积之和

• 这就相当于对于每种分段方式，将其方案数乘上段长之积计入答案

• 令 \(g_{i,j}\) 表示 \(c_i\) 个数分成 \(j\) 段的乘积之和

• 则 \(f_{i,j}\) 可以表示为所有分段方案的价值，乘上这种分段方案强制接上部分分裂点使得总段数变为 \(j\) 的方案数，强制接上的分裂点个数为奇数则扣掉，否则加入

• 可以得到：

• \[f_{i,j}=\sum_{k\ge j}(-1)^{k-j}\binom{k-1}{j-1}g_{i,k} \]

• 在知道 \(g_i\) 的情况下，可以用 FFT 在 \(O(c_i\log c_i)\) 时间内求出 \(f_i\)

• 回到前面如何求 \(g_i\)

• 也就是给定 \(n\) 和 \(m\)，求把 \(n\) 拆分成 \(m\) 个正整数的乘积之和

• 可以视为 \(n+m\) 拆分成 \(m\) 个正整数，所有数 \(-1\) 后的乘积之和

• 也就是在 \(n+m-1\) 个空隙中选 \(m-1\) 个位置插板，把 \(n+m\) 个元素分成 \(m\) 段

• 一段的长度 \(-1\) 的组合意义可以视为在这段中随便选一个空位插板的方案数

• 换句话说，可以看成 \(n+m-1\) 个空隙中选 \(2m-1\) 个位置插板，编号为偶数的板用来切分 \(m\) 段，编号为奇数的板用来在一段中选出一个空隙

• 故把 \(n\) 拆分成 \(m\) 个正整数的乘积之和为 \(\binom{n+m-1}{2m-1}\)！

• 依然对 \(f\) 进行分治 NTT，\(O(m\log m\log n)\)，\(m=\sum_{i=1}^nc_i\)

三、链拓展到环

• 考虑破环为链，设环首位的颜色为 \(x\)

• 这相当于选择恰好一种颜色 \(x\)，强制链的第一段为颜色 \(x\)，并且链的最后一段不能为颜色 \(x\)

• 并且如果要由链构成环，则链的第一段（长度为 \(l\)）中可以选出左边的任意 \([0,l)\) 个元素循环位移到最右边

• 也就是链的第一段有一个为段长的乘积贡献

• 如果不需要限制链的最后一段不为颜色 \(x\) 也不需要考虑第一段的乘积贡献，则令一个新的生成函数 \(G_i(x)\)，与 \(F_i(x)\) 相比只有 \(\frac 1{i!}\) 变成了 \(\frac 1{(i-1)!}\)

• 同时分治 NTT 时计算的是在 \([1,n]\) 中选出恰好一个 \(i\) 为 \(G\) 其他都为 \(F\) 的乘积之和，可以在分治结构上记录两个多项式分别表示是否选出一个 \(G\) 的多项式乘积

• 并且算答案的时候系数乘上的数也是 \((i-1)!\) 而不是 \(i!\)

• 而如果考虑第一段的乘积贡献，\(G_i(x)\) 的系数需要表示的就变成了形如这样的形式：

• \(n\) 拆分成 \(m\) 个数之和的所有方案下，第一个数的二次方乘上其他数的乘积之和

• 和前面一样，把 \(x^2\) 拆成 \(2\binom x2+x\)，\(\binom x2\) 可以给一个组合意义表示第一段中插入 \(2\) 个板的方案数

• 可以得出这个问题的答案为 \(2\binom{n+m-1}{2m}+\binom{n+m-1}{2m-1}\)

• 如果要限制链的最后一段不为颜色 \(x\)，也就是我们需要扣掉首尾都为 \(x\) 的方案数，对应的 \((i-1)!\) 再改成 \((i-2)!\) 再进行一遍分治 NTT 即可

• 为了只做一遍分治 NTT，我们可以把只限制第一段和限制第一段和最后一段的方案放到一起算

• 具体地把 \([x^j]G_i(x)\) 设为：

• \[\frac{w_{i,j}-w_{i,j+1}}{(j-1)!} \]

• 其中 \(w_{i,j}\) 表示把 \(c_i\) 分成若干段，强制合并一些段使得段数变成 \(j\) 的所有段长（第一段平方，其他段不平方，算的是强制合并之前）之积，再乘上 \(-1\) 的强制合并的次数次幂的和

• 可以得出如果设

• \[H(x)=\sum_{i=1}^nG_i(x)\prod_{j\ne i}F_i(x) \]

• 则答案为：

• \[\sum_i(i-1)![x^i]H(x) \]

• 注意这时 \((i-1)![x^i]H(x)\) 的含义是有 \(i\) 段的总方案数（不管最后一段）减去有 \(i+1\) 段的不合法方案数（最后一段和第一段同色）

• \(O(m\log m\log n)\)

Code

#include <bits/stdc++.h>

template <class T>
inline void read(T &res)
{
	res = 0; bool bo = 0; char c;
	while (((c = getchar()) < '0' || c > '9') && c != '-');
	if (c == '-') bo = 1; else res = c - 48;
	while ((c = getchar()) >= '0' && c <= '9')
		res = (res << 3) + (res << 1) + (c - 48);
	if (bo) res = ~res + 1;
}

const int N = 2e5 + 5, M = 1e6 + 5, djq = 998244353;

int n, m, c[N], fac[M], inv[M], f[M], g[M], h[M], ff, tot, rev[M], yg[M], ans;

int C(int n, int m)
{
	return 1ll * fac[n] * inv[m] % djq * inv[n - m] % djq;
}

int qpow(int a, int b)
{
	int res = 1;
	while (b)
	{
		if (b & 1) res = 1ll * res * a % djq;
		a = 1ll * a * a % djq;
		b >>= 1;
	}
	return res;
}

inline void add(int &a, const int &b) {if ((a += b) >= djq) a -= djq;}

inline void sub(int &a, const int &b) {if ((a -= b) < 0) a += djq;}

void FFT(int n, int *a, int op)
{
	for (int i = 0; i < n; i++) if (i < rev[i]) std::swap(a[i], a[rev[i]]);
	int g = qpow(3, (djq - 1) / n * ((n + op) % n));
	yg[0] = 1;
	for (int i = 1; i < n; i++) yg[i] = 1ll * yg[i - 1] * g % djq;
	for (int k = 1, t = n >> 1; k < n; k <<= 1, t >>= 1)
		for (int i = 0; i < n; i += k << 1)
		{
			int w = 1;
			for (int j = 0, x = 0; j < k; j++, x += t)
			{
				int u = a[i + j], v = 1ll * yg[x] * a[i + j + k] % djq;
				add(a[i + j] = u, v); sub(a[i + j + k] = u, v);
			}
		}
	if (op == -1)
	{
		int gg = qpow(n, djq - 2);
		for (int i = 0; i < n; i++) a[i] = 1ll * a[i] * gg % djq;
	}
}

struct poly
{
	std::vector<int> f, g;
} a[N];

std::vector<int> plus(std::vector<int> a, std::vector<int> b)
{
	std::vector<int> res;
	for (int i = 0; i < a.size(); i++)
		res.push_back(a[i]), add(res[i], b[i]);
	return res;
}

std::vector<int> mul(std::vector<int> a, std::vector<int> b)
{
	ff = 1; tot = 0; int n = a.size(), m = b.size();
	while (ff < n + m - 1) ff <<= 1, tot++;
	for (int i = 0; i < ff; i++)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << tot - 1);
	for (int i = 0; i < ff; i++)
		f[i] = i < n ? a[i] : 0, g[i] = i < m ? b[i] : 0;
	FFT(ff, f, 1); FFT(ff, g, 1);
	for (int i = 0; i < ff; i++) f[i] = 1ll * f[i] * g[i] % djq;
	FFT(ff, f, -1);
	std::vector<int> res;
	for (int i = 0; i < n + m - 1; i++) res.push_back(f[i]);
	return res;
}

poly djqtxdy(int l, int r)
{
	if (l == r) return a[l];
	int mid = l + r >> 1;
	poly ls = djqtxdy(l, mid), rs = djqtxdy(mid + 1, r);
	return (poly) {mul(ls.f, rs.f), plus(mul(ls.f, rs.g), mul(ls.g, rs.f))};
}

int main()
{
	read(n);
	for (int i = 1; i <= n; i++) read(c[i]), m += c[i];
	fac[0] = inv[0] = inv[1] = 1;
	for (int i = 1; i <= 1000000; i++) fac[i] = 1ll * fac[i - 1] * i % djq;
	for (int i = 2; i <= 1000000; i++)
		inv[i] = 1ll * (djq - djq / i) * inv[djq % i] % djq;
	for (int i = 2; i <= 1000000; i++) inv[i] = 1ll * inv[i] * inv[i - 1] % djq;
	for (int I = 1; I <= n; I++)
	{
		for (int i = 0; i < c[I]; i++)
		{
			f[c[I] - 1 - i] = 1ll * fac[i] * C(c[I] + i, i * 2 + 1) % djq;
			h[c[I] - 1 - i] = (2ll * fac[i] * C(c[I] + i, i * 2 + 2)
				+ f[c[I] - 1 - i]) % djq;
			if (g[i] = inv[i], i & 1) g[i] = djq - g[i];
		}
		ff = 1; tot = 0;
		while (ff <= (c[I] - 1 << 1)) ff <<= 1, tot++;
		for (int i = 0; i < ff; i++)
			rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << tot - 1);
		for (int i = c[I]; i < ff; i++) f[i] = g[i] = h[i] = 0;
		FFT(ff, f, 1); FFT(ff, g, 1); FFT(ff, h, 1);
		for (int i = 0; i < ff; i++) f[i] = 1ll * f[i] * g[i] % djq,
			h[i] = 1ll * h[i] * g[i] % djq;
		FFT(ff, f, -1); FFT(ff, h, -1);
		for (int i = 0, j = c[I] - 1; i < j; i++, j--)
			std::swap(f[i], f[j]), std::swap(h[i], h[j]);
		for (int i = c[I]; i >= 1; i--) f[i] = 1ll * f[i - 1] * inv[i - 1] % djq,
			h[i] = 1ll * h[i - 1] * inv[i - 1] % djq;
		a[I].f.push_back(0); a[I].g.push_back(0);
		for (int i = 1; i <= c[I]; i++) a[I].f.push_back(1ll * f[i] * inv[i] % djq),
			a[I].g.push_back(1ll * h[i] * inv[i - 1] % djq),
				a[I].g[i - 1] = (a[I].g[i - 1] - 1ll * h[i] *
					(i > 1 ? inv[i - 2] : 0) % djq + djq) % djq;
	}
	poly res = djqtxdy(1, n);
	for (int i = 2; i <= m; i++) ans = (1ll * fac[i - 1] * res.g[i] + ans) % djq;
	return std::cout << ans << std::endl, 0;
}