「雅礼集训 2017 Day11」PATH

发表于
https://loj.ac/problem/6051

题目链接

「雅礼集训 2017 Day11」PATH

做法

$$
概率 = \frac{合法方案数}{总方案数}\\
m = \sum_{i = 1}^{n} a_i\\
总方案数 = \frac{m!}{\prod a_i!}\\
合法方案数 = \frac{m!}{\prod hook(i, j)} = \prod_{1 \leq k \leq i \leq n} \prod_{a_{i + 1} < j \leq a_i} (a_k - k + i - j + 1)\\
= \prod_{1 \leq j \leq i \leq n} \frac{(a_j - a_{i + 1} + i - j)!}{(a_j - a_i + i - j)!}\\
= \frac{\prod_{i = 1}^{n}(n - i + a_i)}{\prod_{1 \leq i < j \leq n}((a_i - i) - (a_j - j))}
$$

$$
概率 = (\prod_{i = 1}^n \frac{a_i!}{(n + a_i - i)!})(\prod_{1 \leq i < j \leq n} ((a_i - i) - (a_j - j)))
$$
我们可以快速处理左边的式子，对于右边的式子，令 $ b_i = a_i - i $ 。

$$
\prod_{1 \leq i \leq j \leq n} ((a_i - i) - (a_j - j)) = \prod_{1 \leq i < j \leq n} (b_i - b_j)
$$

将括号展开发现 $ b_i - b_j > 0 (i < j) $ ，所以当 $ i \geq j $ 时贡献记为负数，可以忽略。

令 $ \prod_{1 \leq i < j \leq n} (b_i - b_j) = \sum x^{f(x)} $ ，显然可以通过卷积求出 $ \sum x^{f(x)} $ 。

注意根据费马小定理 $ f(x) $ 应由 $ mod - 1 $ 取模，所以不能写 $ ntt $ ，由于答案不会爆 $ long~long $ ，所以可以写 $ fft $ 。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
#include <bits/stdc++.h>
#define pb push_back
#define mp make_pair
#define fst first
#define snd second
using namespace std;
typedef long long ll;
typedef long double ld;
const ld pi = acos(-1.0);
const int mod = 1004535809;
const int N = 1000010;
int n, a[N], b[N], ans = 1;
int fac[N], ifac[N];
vector<ll> f, g;

inline int add(const int &x, const int &y) { return x + y < mod ? x + y : x + y - mod; }
inline int sub(const int &x, const int &y) { return x - y < 0 ? x - y + mod : x - y; }
inline int mul(const int &x, const int &y) { return (int)((ll)x * y % mod); }
inline int ksm(int x, ll y = mod - 2) {
	int ss = 1; for(; y; y /= 2, x = mul(x, x)) if(y % 2) ss = mul(ss, x);
	return ss;
}
int C(int x, int y) {
	if(x < 0 || y < 0 || x < y) return 0;
	return mul(fac[x], mul(ifac[y], ifac[x - y]));
}

namespace Poly {
	const int M = 32768;
	inline int Get(int x) { int ss = 1; for(; ss <= x; ss <<= 1); return ss; }
	struct cp { ld x, y; cp(ld X = 0, ld Y = 0) : x(X), y(Y) {} };
	inline cp operator+(cp a, cp b) { return cp(a.x + b.x, a.y + b.y); }
	inline cp operator-(cp a, cp b) { return cp(a.x - b.x, a.y - b.y); }
	inline cp operator*(cp a, cp b) {
		return cp(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
	}
	void dft(vector<cp> &A, int lmt, int opt) {
		A.resize(lmt);
		for(int i = 0, j = 0; i < lmt; i++) {
			if(i > j) swap(A[i], A[j]);
			for(int k = lmt >> 1; (j ^= k) < k; k >>= 1);
		}
		vector<cp> w(lmt >> 1);
		for(int mid = 1; mid < lmt; mid <<= 1) {
			w[0] = cp(1, 0);
			cp wn = cp(cos(opt * pi / mid), sin(opt * pi / mid));
			for(int j = 1; j < mid; j++) w[j] = w[j - 1] * wn;
			for(int R = mid << 1, j = 0; j < lmt; j += R)
				for(int k = 0; k < mid; k++) {
					cp x = A[j + k], y = A[j + mid + k] * w[k];
					A[j + k] = x + y, A[j + mid + k] = x - y;
				}
		}
		if(opt == -1) for(int i = 0; i < lmt; i++) A[i].x /= lmt, A[i].x += 0.5;
	}
	vector<ll> Mul(const vector<ll> &A, const vector<ll> &B) {
		int lmt = Get(A.size() + B.size() - 2);
		vector<cp> a, b; vector<ll> res;
		for(int i = 0; i < A.size(); i++) a.pb(cp(A[i], 0));
		for(int i = 0; i < B.size(); i++) b.pb(cp(B[i], 0));
		dft(a, lmt, 1), dft(b, lmt, 1);
		for(int i = 0; i < lmt; i++) a[i] = a[i] * b[i];
		dft(a, lmt, -1);
		for(int i = 0; i < A.size() + B.size(); i++) res.pb((ll)a[i].x);
		return res;
	}
}

int main() {
	scanf("%d", &n);
	for(int i = 1; i <= n; i++) scanf("%d", &a[i]), b[i] = a[i] + n - i;
	fac[1] = fac[0] = ifac[1] = ifac[0] = 1;
	for(int i = 2; i < N; i++) fac[i] = mul(fac[i - 1], i);
	for(int i = 2; i < N; i++) ifac[i] = mul(mod - mod / i, ifac[mod % i]);
	for(int i = 2; i < N; i++) ifac[i] = mul(ifac[i - 1], ifac[i]);
	for(int i = 1; i <= n; i++) ans = mul(ans, mul(fac[a[i]], ifac[b[i]]));
	if(n <= 3000) {
		for(int i = 1; i < n; i++)
			for(int j = i + 1; j <= n; j++)
				ans = mul(ans, a[i] - a[j] + j - i);
		printf("%d\n", ans); return 0;
	}
	f.resize(N), g.resize(N);
	for(int i = 1; i <= n; i++) ++f[a[i] - i + n], ++g[i - a[i] + a[1]];
	for(; !f.back(); f.pop_back()); for(; !g.back(); g.pop_back());
	f = Poly::Mul(f, g), f.resize(N);
	for(int i = a[1] + n + 1; i <= (a[1] + n) * 2; i++)
		ans = mul(ans, ksm(i - a[1] - n, f[i]));
	printf("%d\n", ans);
	return 0;
}