【组合数学】【P5216】DLS采花

Description

DLS 有 (N) 个花田，每个花田里有 (a_i) 朵花。

DLS 喜欢稀奇古怪的花田，他希望重新排列花田，然后去采花。
但 DLS 采花又有一个癖好：他会从左往右采花。
若当前采到第 (i) 个花田，在之前有一个花田的花的数量，是第 (i) 个花田的花的数量的因子的话，那么 DLS 不会采这个花田的花。

现在，DLS 想知道对于所有排列花田的方案，他能够采到的花的数量的和是多少。
由于答案会比较大，请对 (998244353) 取模。

Input

第一行是一个整数 (N)

下面一行 (N) 个整数代表这个序列

Output

输出答案对 (998244353) 取膜的结果

Hint

(for~50~percents:N~leq~9)

(for~80~percent:N~leq~1000)

(forall:1~leq~N~leq~10^5,1~leq~10^5)

Solution

感觉所有的题解都没有写部分分和式子是怎么用数学方法推出来的啊qwq

考虑前50分，直接枚举全排列检验，时间复杂度 (O(n!~ imes~n^2))，期望得分50pts。

考虑前80的数据，我们发现问题等价于求(sum_{i = 1}^n a_i ~ imes~ ext{序列中所有}a_i ext{的因数都在}a_i ext{的后面的方案数})。问题转化为如何求出方案数。

对于每个 (a_i)，我们考虑枚举他在产生贡献方案中的位置。设 (a_i) 除去他自身共有 (k) 个因数在序列中。

当 (a_i) 放在倒数第 (k + 1) 位置上时，后面 (k) 个因数随便放，前面的所有数字随便放，所以方案数为 (k!~ imes~(n - k - 1)!)

当 (a_i) 放在倒数第 (k + 2) 位置上时，后面 (k) 个因数和随便一个非因数随便放，前面剩下的数字随便放。由于非因数是随便选的，所以方案数为 ((k + 1)!~ imes~(n - k - 2)!~ imes~C_{n - k - 1}^{1})

当 (a_i) 放在倒数第 (k + 3) 位置上时，后面 (k) 个因数和随便两个非因数随便放，前面剩下的数字随便放。由于非因数是随便选的，所以方案数为 ((k + 2)!~ imes~(n - k - 3)!~ imes~C_{n - k - 1}^{2})

一位读者砸烂了复读机停止了复读

依此做数学归纳，可得答案即为

[sum_{i = k}^{n - 1} i!~ imes~(n - i + 1)!~ imes~C_{n - k - 1}^{i - k} ]

相当于枚举 (a_i) 在倒数第 (i + 1) 位时的答案。

发现这个式子sigma后面可以在预处理阶乘和逆元后 (O(1)) 计算，整个sigma可以 (O(n)) 计算，一共算 (n) 次，于是总复杂度 (O(n^2))，期望得分80pts。

考虑全部的数据，我们发现上面先枚举 (a_i) 的位置再枚举因数怎么放的多余的，可以直接把他们放在一起枚举：我们在一个序列中选择 ((k + 1)) 个位置，其中第一个位置放 (a_i)，其他位置放因数，剩下的位置随便放，于是答案即为

[C_{n}^{k + 1}~ imes~k!~ imes~(n - k - 1)! ]

这个式子在预处理后可以 (O(1)) 计算，于是算上求因数后总复杂度 (O(n~sqrt a))，期望得分100pts。

Code

#include <cstdio>
#ifdef ONLINE_JUDGE
#define freopen(a, b, c)
#define printtime()
#else
#include <ctime>
#define printtime() printf("Times used  = %ld ms
", clock())
#endif
#define ci const int
#define cl const long long

typedef long long int ll;

namespace IPT {
	const int L = 1000000;
	char buf[L], *front=buf, *end=buf;	
	char GetChar() {
		if (front == end) {
			end = buf + fread(front = buf, 1, L, stdin);
			if (front == end) return -1;
		}
		return *(front++);
	}
}

template <typename T>
inline void qr(T &x) {
	char ch = IPT::GetChar(), lst = ' ';
	while ((ch > '9') || (ch < '0')) lst = ch, ch=IPT::GetChar();
	while ((ch >= '0') && (ch <= '9')) x = (x << 1) + (x << 3) + (ch ^ 48), ch = IPT::GetChar();
	if (lst == '-') x = -x;
}

namespace OPT {
	char buf[120];
}

template <typename T>
inline void qw(T x, const char aft, const bool pt) {
	if (x < 0) {x = -x, putchar('-');}
	int top=0;
	do {OPT::buf[++top] = static_cast<char>(x % 10 + '0');} while (x /= 10);
	while (top) putchar(OPT::buf[top--]);
	if (pt) putchar(aft);
}

const int maxn = 100010;
const int MOD = 998244353;

int n;
int MU[maxn], bk[maxn];
ll ans;
ll fact[maxn], inv[maxn], finv[maxn];

ll C(int, int);

int main() {
	freopen("1.in", "r", stdin);
	qr(n);
	for (int i = 1; i <= n; ++i) {
		qr(MU[i]); ++bk[MU[i]];
	}
	fact[0] = finv[0] = inv[1] = 1;
	for (int i = 2; i <= n; ++i) inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
	for (int i = 1; i <= n; ++i) fact[i] = fact[i - 1] * i % MOD, finv[i] = finv[i - 1] * inv[i] % MOD;
	for (int i = 1; i <= n; ++i) {
		int _cnt = 0;
		int j;
		for (j = 1; (j * j) < MU[i]; ++j) if (!(MU[i] % j)) _cnt += bk[j] + bk[MU[i] / j];
		if ((j * j) == MU[i]) _cnt += bk[j];
		--_cnt;
		ans = (ans + C(n, _cnt + 1) * fact[_cnt] % MOD * fact[n - _cnt - 1] % MOD * MU[i]) % MOD;
	}
	qw((ans + MOD) % MOD, '
', true);
	printtime();
}

ll C(int x, int y) {
	return fact[x] * finv[y] % MOD * finv[x - y] % MOD;
}