HAOI2018 简要题解

这套题是 dy, wearry 出的。学长好强啊，可惜都 (wc) 退役了。。

话说 wearry 真的是一个计数神仙。。就没看到他计不出来的题。。。每次考他模拟赛总有一两道毒瘤计数TAT

上午的官方题解可以看 dy0607 的博客，写的挺详细的。

「HAOI2018」奇怪的背包

题意

小C非常擅长背包问题，他有一个奇怪的背包，这个背包有一个参数 (P) ，当他向这个背包内放入若干个物品后，背包的重量是物品总体积对 (P) 取模后的结果．

现在小C有 (n) 种体积不同的物品，第 (i) 种占用体积为 (V_i) ，每种物品都有无限个．他会进行 (q) 次询问，每次询问给出重量 (w_i) ，你需要回答有多少种放入物品的方案，能将一个初始为空的背包的重量变为 (w_i) ．注意，两种方案被认为是不同的，当且仅当放入物品的种类不同，而与每种物品放入的个数无关．不难发现总的方案数为 (2^n) .

由于答案可能很大，你只需要输出答案对 (10^9 + 7) 取模的结果．

(1 le n, q le 10^6, 3 le P le 10^9, 0 < V_i, w_i < P)

题解

题目就是求，有多少种选择物品的方式，使得下面的同余方程有解，（假设选择的物品的体积为 (a_1, a_2, dots, a_k) ）

[sum_{i = 1}^k x_ia_i equiv w_j pmod{P}\ Rightarrow sum_{i = 1}^k x_ia_i + x_0P = w_j ]

由于裴蜀定理，有解当且仅当 (gcd(a_1, a_2, …, a_k, P) mid w_j) 。

然后令 (f_{i, x}) 为前 (i) 种物品，选择物品的体积与 (P) 做 (gcd) 结果为 (x) 的方案数。

直接转移是 (mathcal O(n sigma(P) log P)) 的，(sigma(P)) 为 (P) 的约数个数，大概是 (P^{frac 13}) 的级别。

这样显然还不够优秀，不难发现对于 (gcd(a_i, P)) 相同的 (i) 的转移其实是一样的，我们把这些合并一起转移就好啦。

此时复杂度优化到了 (mathcal O((sigma^2(P) + n + q) log P)) 其中 (log P) 为取 (gcd) 的复杂度。

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("2523.in", "r", stdin);
	freopen ("2523.out", "w", stdout);
#endif
}

const int N = 5100, Mod = 1e9 + 7;

int n, q, P, d[N], cnt;

inline int id(int x) {
	return lower_bound(d + 1, d + cnt + 1, x) - d;
}

int dp[N][N], tot[N], Pow[1001000], ans[N];

int main () {

	File();

	n = read(); q = read(); P = read();

	For (i, 1, sqrt(P + .5)) if (!(P % i)) {
		d[++ cnt] = i;
		if (i * i != P) d[++ cnt] = P / i;
	}
	sort(d + 1, d + cnt + 1);

	For (i, 1, n)
		++ tot[id(__gcd(P, read()))];

	Pow[0] = 1;
	For (i, 1, n) Pow[i] = Pow[i - 1] * 2 % Mod;

	dp[0][cnt] = 1;
	For (i, 1, cnt) For (j, 1, cnt) {
		int x = id(__gcd(d[i], d[j]));
		dp[i][x] = (dp[i][x] + 1ll * dp[i - 1][j] * (Pow[tot[i]] - 1)) % Mod;
		dp[i][j] = (dp[i][j] + dp[i - 1][j]) % Mod;
	}

	For (i, 1, cnt) For (j, 1, cnt) 
		if (!(d[i] % d[j])) (ans[i] += dp[cnt][j]) %= Mod;

	while (q --) printf ("%d
", ans[id(__gcd(read(), P))]);

    return 0;

}

「HAOI2018」反色游戏

题意

小C和小G经常在一起研究搏弈论问题，有一天他们想到了这样一个游戏．

有一个 (n) 个点 (m) 条边的无向图，初始时每个节点有一个颜色，要么是黑色，要么是白色．现在他们对于每条边做出一次抉择：要么将这条边连接的两个节点都反色（黑变白，白变黑），要么不作处理．他们想把所有节点都变为白色，他们想知道在 (2^m) 种决策中，有多少种方案能达成这个目标．

小G认为这个问题太水了，于是他还想知道，对于第 (i) 个点，在删去这个点及与它相连的边后，新的答案是多少．

由于答案可能很大，你只需要输出答案对 (10^9 + 7) 取模后的结果．

对于所有数据，有 (1 le T le 5, 1 le n, m le 10^5, 1 le u, v le n) ，没有重边和自环。

题解

一道有意思的 找规律 结论题 qwq

首先我们不管删点操作，思考一下什么时候是无解的。一个显然的必要条件是，一个连通块内必须要有偶数个黑点，而这个条件也是充分的。为什么呢？

我们将黑点任意两两配对，对一对黑点考虑它们之间的任意一条路径，将路径上的边的状态（操作or不操作）取反，立即可以得到一种合法方案。

同时我们发现，对于任意一种偶数个黑点的方案，如果进行上面的配对，都可以唯一对应到一种全为白色的方案。

也就是说，只要有解，解的个数都是相同的。对于一个连通块，假设边数为 (m) ，点数为 (n) ，有 (2^{n − 1}) 种对点染色的方案是有解的且解的个数相同，总方案数为 (2^m)，那么对于一个有解的局面这个连通块的方案数就是 (2^{m − n + 1}) ! 对于整张图，总的方案数就是 (2^{m − n + p}) ，其中 (p) 是连通块个数。

看起来不太好想，规律其实 很好找 。

有一种用线性基理解的方式更加优秀，不难发现每条边对应一个长度为 (n) 的位向量，其中只有两个点为 (1) 。对于一个联通块来说，如果有解就是要异或出 (0) ，然后方案数是 (2^{free}) 。

不难发现自由元 (free) 的个数其实就是 (m - n + 1) ，因为极大线性无关组大小为 (n - 1) 是一定可达的，然后一开始共有 (m) 个向量，减一减就是的了。

然后剩下只需要维护每个点被删去的答案，这个用类似于 Tarjan 求割点的方法去做。具体来说只需要维护一下断开某个点会产生几个联通块，以及每个联通块是否合法即可。

复杂度是 (mathcal O(n + m)) 的。

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("2524.in", "r", stdin);
	freopen ("2524.out", "w", stdout);
#endif
}

const int N = 1e5 + 1e3, Mod = 1e9 + 7;

int n, m, deg[N]; char S[N];

vector<int> G[N]; bool pass[N]; 

int dfn[N], lowlink[N], clk;
int cut[N], sz[N], sub[N], bel[N], now;

#define chk(x) !((x) & 1)

void Tarjan(int u, int fa = 0) {
	dfn[u] = lowlink[u] = ++ clk;
	pass[u] = true; sz[u] = (S[u] ^ 48); 
	bel[u] = now; sub[u] = cut[u] = 0;
	for (int v : G[u]) if (!dfn[v]) {
		Tarjan(v, u);
		chkmin(lowlink[u], lowlink[v]);
		sz[u] += sz[v];
		if (lowlink[v] >= dfn[u]) {
			++ cut[u];
			sub[u] += sz[v];
			pass[u] &= chk(sz[v]);
		}
	} else if (v != fa) chkmin(lowlink[u], dfn[v]);
	if (!fa) -- cut[u];
}

int Pow2[N];

int main () {

	File();

	int cases = read();
	while (cases --) {

		n = read(); m = read();
		For (i, 1, n) G[i].clear(), deg[i] = dfn[i] = 0;

		For (i, 1, m) {
			int u = read(), v = read();
			G[u].push_back(v); G[v].push_back(u);
			++ deg[u]; ++ deg[v];
		}
		scanf ("%s", S + 1);

		int cnt_odd = 0, ans = m - n;
		For (i, 1, n) if (!dfn[i])
			now = i, Tarjan(i), cnt_odd += (sz[i] & 1), ++ ans;

		Pow2[0] = 1;
		For (i, 1, m) Pow2[i] = Pow2[i - 1] * 2 % Mod;

		printf ("%d", cnt_odd ? 0 : Pow2[ans]);
		For (i, 1, n)
			if (!deg[i]) {
				printf (" %d", (cnt_odd - sz[i]) ? 0 : Pow2[ans]);
			} else {
				if (pass[i] && chk(sz[bel[i]] - sub[i] - (S[i] ^ 48)) && !(cnt_odd - (sz[bel[i]] & 1)))
					printf (" %d", Pow2[ans - deg[i] + (cut[i] + 1)]);
				else printf (" 0");
			}
		puts("");

	}

	return 0;

}

「HAOI2018」字串覆盖

题意

小C对字符串颇有研究，他觉得传统的字符串匹配太无聊了，于是他想到了这样一个问题．

对于两个长度为 (n) 的串 (A, B) , 小C每次会给出给出 (4) 个参数 (s, t, l, r) . 令 (A) 从 (s) 到 (t) 的子串(从 (1) 开始标号)为 (T)，令 (B) 从 (l) 到 (r) 的子串为 (P).然后他会进行下面的操作：

如果 (T) 的某个子串与 (P) 相同，我们就可以覆盖 (T) 的这个子串，并获得 (K - i) 的收益，其中 (i) 是初始时 (A) 中(注意不是 (T) 中)这个子串的起始位置，(K)是给定的参数．一个位置不能被覆盖多次．覆盖操作可以进行任意多次，你需要输出获得收益的最大值．

注意每次询问都是独立的，即进行一次询问后，覆盖的位置会复原．

对于所有数据，有 (1 le n, q le 10^5) ，(A, B) 仅由小写英文字母组成，(1 le s le t le n, 1 le l le r le n, n < K le 10^9)．

对于 (n = 10^5) 的测试点，满足 (51 le r - l le 2000) 的询问不超过 (11000) 个，且 (r - l) 在该区间内均匀随机．

题解

其实那个数据范围提示了一下做法。。。就是对于 (r - l) 的长度进行讨论。

对于 (r - l) 较大的点，显然段数比较少，这个地方用数据结构模拟选取方案就行了。一个比较经典的解法是 后缀数组+主席树 。

我们把 (underline{A#B}) 做后缀排序，就是我们可以选取的串是和 (B_{l, r}) 相同的串，这对应了后缀数组上的一段连续区间。然后我们要从这个区间里选取尽量靠前的点，满足和前面选取的区间不重复。这个直接在主席树上二分即可。

那么现在只需要解决 (r - l) 比较小的情况。我们对于每种长度单独考虑即可。就是从前往后依次选很多段，每个相邻两端连边的话就构成了一个树上结构，利用树上倍增统计答案即可。

利用均值法在 (sqrt n) 分开讨论复杂度是最优的，为 (mathcal O(n sqrt n log n)) 。

但是这样应该会跑的比较慢，dy 为了这个显著快于暴力，限制了长度的频率，我们在 (50) 分开就行了。

好像 孔爷 在知乎上回答可以做到 (mathcal O(n sqrt n)) 。

具体来说继续利用均值法，前面用值域分块，后面用轻重链剖分贪心。

代码

实在比较难写。。抄了一些 dy 的写法， dyy 码力还是强啊。

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;
using ll = long long;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("2525.in", "r", stdin);
	freopen ("2525.out", "w", stdout);
#endif
}

const int N = 2e5 + 1e3, Gap = 50;

int Log2[N], n, K, q; char S[N]; 

int m, sa[N], rk[N], tmp[N], pre[N];

void Radix_Sort() {
	For (i, 1, m) pre[i] = 0;
	For (i, 1, n) ++ pre[rk[i]];
	For (i, 1, m) pre[i] += pre[i - 1];
	Fordown (i, n, 1) sa[pre[rk[tmp[i]]] --] = tmp[i];
}
void Build_SA() {
	For (i, 1, n) rk[tmp[i] = i] = S[i];
	m = 255; Radix_Sort();
	for (int k = 1, p = 0; k <= n; k <<= 1, p = 0) {
		For (i, n - k + 1, n) tmp[++ p] = i;
		For (i, 1, n) if (sa[i] > k) tmp[++ p] = sa[i] - k;
		Radix_Sort(); swap(rk, tmp);
		rk[sa[1]] = m = 1;
		For (i, 2, n)
			rk[sa[i]] = (tmp[sa[i]] == tmp[sa[i - 1]] && tmp[sa[i] + k] == tmp[sa[i - 1] + k]) ? m : ++ m;
		if (m >= n) break;
	}
}

int height[N];
void Get_Height() {
	for (int i = 1, j, k = 0; i <= n; ++ i) {
		if (k) -- k; j = sa[rk[i] + 1];
		if (!j) continue;
		while (S[i + k] == S[j + k]) ++ k;
		height[rk[i]] = k;
	}
}

template<int Maxn>
struct President_Segment_Tree {

	int ls[Maxn], rs[Maxn], tot[Maxn], Node;

#define mid ((l + r) >> 1)

	void Update(int &o, int pre, int l, int r, int up) {
		tot[o = ++ Node] = tot[pre] + 1; ls[o] = ls[pre]; rs[o] = rs[pre];
		if (l == r) return;
		if (up <= mid) Update(ls[o], ls[pre], l, mid, up);
		else Update(rs[o], rs[pre], mid + 1, r, up);
	}

	int Find(int x, int y, int l, int r, int qp) {
		if (!(tot[y] - tot[x])) return 0;
		if (l == r) return l;
		int res = qp <= mid ? Find(ls[x], ls[y], l, mid, qp) : 0;
		return res ? res : Find(rs[x], rs[y], mid + 1, r, qp);
	}

};

President_Segment_Tree<N * 20> T; int rt[N];

int minv[N][22];

struct Query {
	int s, t, l, len, id, pos;
} Q[N];

int fa[N][22]; ll sum[N][22], ans[N];

inline void Work(int len, int &it) {
	static int pos[N]; 
	if (Q[it].len != len) return;
	for (int l = 1, r = l; l <= n * 2 + 1; r = l = r + 1) {
		int cnt = 0;
		while (height[r] >= len) ++ r;
		For (j, l, r) if (sa[j] <= n) pos[++ cnt] = sa[j];

		if (Q[it].pos < l || Q[it].pos > r) continue;

		sort(pos + 1, pos + cnt + 1);
		int lim = Log2[min(cnt + 1, n / len)], k = 1;
		pos[cnt + 1] = 2 * n + 1;

		For (i, 1, cnt) {
			while (pos[k] - pos[i] < len) ++ k;
			fa[i][0] = k, sum[i][0] = K - pos[i];
		}

		For (j, 1, lim) For (i, 1, cnt) {
			fa[i][j] = fa[fa[i][j - 1]][j - 1];
			sum[i][j] = sum[i][j - 1] + sum[fa[i][j - 1]][j - 1];
		}

		for (; Q[it].len == len && Q[it].pos <= r; ++ it) {
			int s = Q[it].s, t = Q[it].t, x = lower_bound(pos + 1, pos + cnt + 2, s) - pos;
			if (pos[x] <= t) {
				int u = x; ll res = 0;
				Fordown (i, lim, 0) if (pos[fa[u][i]] && pos[fa[u][i]] <= t)
					res += sum[u][i], u = fa[u][i];
				ans[Q[it].id] = res + sum[u][0];
			}
		}

		For (i, 1, cnt + 1) For (j, 0, lim) fa[i][j] = sum[i][j] = 0;
	}
}

int main () {

	File();

	n = read(); K = read();

	scanf ("%s", S + 1); S[n + 1] = '#';
	scanf ("%s", S + n + 2);

	n = n * 2 + 1; Build_SA(); Get_Height(); 

	For (i, 1, n) minv[i][0] = height[i];
	For (i, 2, n) Log2[i] = Log2[i >> 1] + 1;
	For (j, 1, Log2[n])
		for (int i = 1; i + (1 << j) <= n; ++ i)
			minv[i][j] = min(minv[i][j - 1], minv[i + (1 << (j - 1))][j - 1]);

	n >>= 1;

	q = read();
	For (i, 1, q) {
		int s = read(), t = read(), l = read(), r = read();
		Q[i] = (Query) {s, t - (r - l), l, r - l + 1, i, rk[l + n + 1]};
	}
	sort(Q + 1, Q + q + 1, [&](Query lhs, Query rhs) { return lhs.len != rhs.len ? lhs.len < rhs.len : lhs.pos < rhs.pos; } );

	int it = 1;
	For (i, 1, min(n, Gap)) Work(i, it);
	For (i, 1, 2 * n + 1) {
		if (sa[i] <= n) T.Update(rt[i], rt[i - 1], 1, n, sa[i]);
		else rt[i] = rt[i - 1];
	}

	for (; it <= q; ++ it) {
		Query A = Q[it];
		int x = A.pos, l = x, r = x;
		Fordown (i, Log2[n], 0) 
			if (l - (1 << i) > 0 && minv[l - (1 << i)][i] >= A.len) l -= 1 << i;
		Fordown (i, Log2[n], 0)
			if (r + (1 << i) <= 2 * n + 1 && minv[r][i] >= A.len) r += 1 << i;
		for(int pos = A.s, cur; pos <= n; pos = cur + A.len) {
			cur = T.Find(rt[l - 1], rt[r], 1, n, pos);
			if (!cur || cur > A.t) break; 
			ans[A.id] += K - cur;
		}
	}

	For (i, 1, q) printf ("%lld
", ans[i]);

	return 0;

}

「HAOI2018」苹果树

题意

小C在自己家的花园里种了一棵苹果树，树上每个结点都有恰好两个分支。经过细心的观察，小C发现每一天这棵树都会生长出一个新的结点。

第一天的时候, 果树会长出一个根结点，以后每一天，果树会随机选择一个当前树中没有长出过结点的分支，
然后在这个分支上长出一个新结点，新结点与分支所属的结点之间连接上一条边。

小C定义一棵果树的不便度为树上两两结点之间的距离之和，两个结点之间的距离定义为从一个点走到另一个点的路径经过的边数。

现在他非常好奇，如果 (N) 天之后小G来他家摘苹果，这个不便度的期望 (E) 是多少。但是小C讨厌分数,，所以他只想知道 (E imes N!) 对 (P) 取模的结果，可以证明这是一个整数。

(N le 2000, P le 10^9 + 7)

题解

还是太菜了，这道题还是记不清。

好像有千奇百怪的做法，最简单的还是考虑边对于点对的贡献。

你考虑 (i o fa_i) 这条边被多少个点对穿过，考虑枚举一下 (i) 子树的大小 (j) ，那么会被 (j imes (n - j)) 对点穿过。

我们只需要统计 (i) 子树大小为 (j) 的方案数，先考虑算子树内的方案数，也就是从 (n - i) 个点里选 (j - 1) 个点乘上子树内的合法结构数，也就是 (displaystyle {n - i choose j - 1} imes j!) 。

然后考虑一下子树外的点，(1 sim i) 的排列方式有 (i!) 种，剩下的 (n - i - j + 1) 个点需要逐个插上这棵树，但是不能插到 (i) 号点的子树内，这一部分贡献其实就是 (displaystyle i! imes prod_{k = 1}^{n - i - j + 1} (i + k - 2)) 。

那么最后的答案就是

[ans = sum_{i = 1}^{n} sum_{j = 1}^{n - i + 1} (j imes (n - j)) ({n - i choose j - 1} j!)(i! imes prod_{k = 1}^{n - i - j + 1}(i + k - 2)) ]

(P) 可能没有逆元有点烦，但考虑预处理 mulfac[i][j] 为 (prod_{k = 0}^{j} (i + k)) ，就可以做到 (mathcal O(n^2)) 。

其实有一种更好的方式，把前面那个式子化简就变成

[sumlimits_{i=1}^nsumlimits_{j=1}^{n-i+1}j! (n-j)! {n - i choose j - 1} i(i-1)j ]

就直接能做到 (mathcal O(n^2)) 啦。

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("2526.in", "r", stdin);
	freopen ("2526.out", "w", stdout);
#endif
}

const int N = 2010;

int n, f[N], fac[N], Mod;

int Comb[N][N], mulfac[N][N];

int main () {

	File();

	n = read(); Mod = read();

	fac[0] = 1;
	For (i, 1, n) fac[i] = 1ll * fac[i - 1] * i % Mod;

	For (i, 0, n) {
		Comb[i][0] = 1;
		For (j, 1, i)
			Comb[i][j] = (Comb[i - 1][j] + Comb[i - 1][j - 1]) % Mod;
	}
	For (i, 1, n) {
		mulfac[i][0] = 1;
		For (j, 1, n)
			mulfac[i][j] = 1ll * mulfac[i][j - 1] * (i + j - 1) % Mod;
	}

	int ans = 0;
	For (i, 2, n) For (j, 1, n - i + 1) {
		int coef = mulfac[i - 1][n - (i + j) + 1], tmp = 1;
		ans = (ans + 1ll * fac[j] * Comb[n - i][j - 1] % Mod * (j * (n - j)) % Mod * coef % Mod * fac[i]) % Mod;
	}
	printf ("%d
", ans);

	return 0;

}

「HAOI2018」染色

题意

为了报答小 C 的苹果, 小 G 打算送给热爱美术的小 C 一块画布,
这块画布可以抽象为一个长度为 (N) 的序列, 每个位置都可以被染成 (M) 种颜色中的某一种。

然而小 C 只关心序列的 (N) 个位置中出现次数恰好为 (S) 的颜色种数，如果恰好出现了 (S) 次的颜色有 (K) 种, 则小C会产生 (W_k) 的愉悦度。

小 C 希望知道对于所有可能的染色方案, 他能获得的愉悦度的和对 (1004535809) 取模的结果是多少。

(N le 10^7, M le 10^5, S le 150, 0 le W_i < 1004535809)

题解

看到 恰好 不难想到算 至少 。

我们令 (f_i) 为至少有 (i) 种颜色满足恰好出现 (S) 次的方案数。

颜色选法有 (displaystyle {m choose i}) 种，位置选法有 (displaystyle {n choose iS}) 种，选了的位置内部的排列就是一个可重集元素排列数有 (displaystyle frac{(iS)!}{(S!)^i}) 种，然后其他位置我们随意排除了这 (i) 种颜色的数有 ((m - i)^{n - iS}) （特别的，我们可以令 (0^0 = 1) 虽然好像在数学上没有意义。。）。

那么我们现在算的 (f_i) 就是

[f_i = {m choose i} {n choose iS} frac{(iS)!}{(S!)^i} (m - i)^{n - iS} ]

但是这样算了之后，不难发现 (f_i) 算的答案会有点小问题，就是对于一个排列我们会重复计算很多次。

我们可以考虑 容斥 ，也就是考虑 (f_j(i < j)) 在 (f_i) 中算进去了几次。（假设当前的 (f_j) 取值刚好不重不漏）

不难发现对于每 (i) 个恰好出现 (S) 次的都会被算进去一次，那么总共会被算 (displaystyle {j choose i}) 次。

(displaystyle i = min(m, lfloor frac{n}{S} floor)) 算的 (f_i) 一开始的取值一定正确，因为剩下的位置不够再出现一个恰好 (S) 次的颜色。

那么我们就有一个 (mathcal O(n^2)) 的容斥啦。

[g_i = f_i - sum_{j > i} g_j {j choose i} ]

此时得到的 (g_i) 就是恰好有 (i) 种颜色出现 (S) 次的方案数！

显然是可以用分治 (FFT) 优化的，但这么经典的卷积形式显然没必要。

考虑移项

[f_i = sum_{j ge i} {j choose i} g_i ]

不难发现这是个经典的二项式反演形式（广义容斥定理）。

[g_i = sum_{j ge i} (-1)^{j - i} {j choose i} f_i ]

这样就直接做一遍卷积就好啦，复杂度是 (mathcal O(m log m)) 的。

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
#define div Div

using namespace std;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline int read() {
	int x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("2527.in", "r", stdin);
	freopen ("2527.out", "w", stdout);
#endif
}

const int N = 1e7 + 1e3, Mod = 1004535809;

inline int fpm(int x, int power) {
	int res = 1;
	for (; power; power >>= 1, x = 1ll * x * x % Mod)
		if (power & 1) res = 1ll * res * x % Mod;
	return res;
}

inline void add(int &x, int y) {
	if ((x += y) >= Mod) x -= Mod;
}

inline void sub(int &x, int y) {
	if ((x -= y) < 0) x += Mod;
}

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

inline int div(int x, int y) {
	return 1ll * x * fpm(y, Mod - 2) % Mod;
}

namespace Poly {

	const int Maxn = 1 << 20, g = 3;

	int powg[Maxn], invpowg[Maxn];

	void NTT_Init() {
		for (int i = 2; i < Maxn; i <<= 1)
			invpowg[i] = fpm(powg[i] = fpm(g, (Mod - 1) / i), Mod - 2);
	}

	int len, rev[Maxn];

	void NTT(int *P, int opt) {
		Rep (i, len) if (i < rev[i]) swap(P[i], P[rev[i]]);
		for (int i = 2, p = 1; i <= len; p = i, i <<= 1) {
			int Wi = opt == 1 ? powg[i] : invpowg[i];
			for (int j = 0; j < len; j += i)
				for (int k = 0, x = 1; k < p; ++ k) {
					int u = P[j + k], v = mul(x, P[j + k + p]);
					P[j + k] = (u + v) % Mod; 
					P[j + k + p] = (u - v + Mod) % Mod; 
					x = mul(x, Wi);
				}
		}
		if (!~opt) {
			int inv = fpm(len, Mod - 2);
			Rep (i, len) P[i] = mul(P[i], inv);
		}
	}

	int A[Maxn], B[Maxn], C[Maxn];
	void Mult(int *a, int *b, int *c, int na, int nb) {
		int nc = na + nb, bit = 0;
		For (i, 0, na) A[i] = a[i];
		For (i, 0, nb) B[i] = b[i];

		for (len = 1; len <= nc; len <<= 1) ++ bit;
		Rep (i, len) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
		NTT(A, 1); NTT(B, 1);
		Rep (i, len) C[i] = mul(A[i], B[i]);
		NTT(C, -1);
		For (i, 0, nc) c[i] = C[i];
	}

}

int fac[N], ifac[N];

void Math_Init(int maxn) {
	fac[0] = ifac[0] = 1;
	For (i, 1, maxn) fac[i] = mul(fac[i - 1], i);
	ifac[maxn] = fpm(fac[maxn], Mod - 2);
	Fordown (i, maxn - 1, 1) ifac[i] = mul(ifac[i + 1], i + 1);
}

inline int comb(int n, int m) {
	if (n < 0 || m < 0 || n < m) return 0;
	return mul(mul(fac[n], ifac[n - m]), ifac[m]);
}

int n, m, s, w[N], f[N], g[N];

int main () {

	File();

	n = read(); m = read(); s = read();
	For (i, 0, m) w[i] = read();
	Math_Init(max(n, m)); 
	Poly :: NTT_Init();

	int lim = min(m, n / s);
	For (i, 0, lim) {
		int coef = mul(comb(n, i * s), comb(m, i));
		coef = mul(coef, div(fac[i * s], fpm(fac[s], i)));
		f[i] = mul(coef, fpm(m - i, n - i * s));
	}

	reverse(f, f + lim + 1);
	For (i, 0, lim)
		f[i] = mul(i & 1 ? Mod - 1 : 1, mul(f[i], fac[lim - i]));
	Poly :: Mult(f, ifac, g, lim, lim);
	For (i, 0, lim)
		g[i] = mul(mul(i & 1 ? Mod - 1 : 1, ifac[lim - i]), g[i]);
	reverse(g, g + lim + 1);

	int ans = 0;
	For (i, 0, lim)
		add(ans, mul(g[i], w[i]));
	printf ("%d
", (ans + Mod) % Mod);

	return 0;

}