LOJ6069 「2017 山东一轮集训 Day4」塔

题目来源：LibreOJ #6069. 「2017 山东一轮集训 Day4」塔

题目大意

我们的主人公 yzh 有 (L) 个妹子，她们站成一排，依次编号为 (1dots L)。

一天晚上，yzh 要光♂顾这些妹子。经过一年的养“精”蓄锐，yzh 现在有 (n) 点体力值，每光顾一个妹子会使他的体力值减少 (1)。因为 yzh 很好色，所以他会光顾 (n) 个妹子（不重复），也就是用完所有体力值。

众所周知，男人太快是不好的，而随着 yzh 体力值变少，他就会越来越快。具体来说，假如他光顾某个妹子 (i) 之前，体力值为 (k)，则这个妹子会产生相应的爽感，并发出音量为 (k) 的娇喘。这会恰好被所有距离当前妹子不超过 (k) 的妹子听到，也就是被编号在 ([max(1, i - k),i - 1]cup[i + 1,min(L,i + k)]) 的这些妹子听到。

为了后宫和谐，yzh 不希望任何妹子听到别人的娇喘。也就是说，对于任意两个妹子，假设它们的坐标分别为 (i_1,i_2)，yzh 光临她们时体力值分别为 (k_1,k_2)，则 (|i_1-i_2|>max(k_1,k_2))。

求 yzh 有多少种光临妹子的方案。两种方案不同，当且仅当 yzh 光临的妹子集合不同，或光临顺序不同。方案数对 (m) 取模。

数据范围 (1leq Lleq 10^9)，(1leq nleq 100)，(1leq mleq 10^9)。

本题题解

（强烈建议大家读完有趣的【题目大意】部分）。

设 yzh 光临每个妹子时的体力值（也就是原题面里每座塔的高度）分别为 (p_1,p_2,dots ,p_n)。设 (s = sum_{i = 2}^{n}max(p_i,p_{i - 1}))，则把整个方案紧密地排在一起，需要 (s+1) 个格子。此时还剩下 (L - s - 1) 个格子。把这些格子，放在 (n+1) 个间隙（含两边）里，可以为空，根据插板法，方案数是 ({L - s - 1 + nchoose n})。

(s) 最大不超过 (sum_{i=1}^{n}2i)，即 (n(n+1))。

考虑先对每个 (sin [0,n(n+1)])，求出 (sum_{i = 2}^{n}max(p_i,p_{i - 1}) = s) 的排列 (p) 的数量。注意，(n(n + 1)) 只是我们粗略估计的 (s) 的一个上界，事实是达不到这个上界的，不过不要紧，不存在这样的 (s) 时算方案数为 (0) 即可。

用 DP 求。设 (dp[i][j][k]) 表示从小到大考虑了 (1dots i) 这些数字（已将它们加入排列），当前 (s) 的值为 (j)，已加入的数字形成了 (k) 个连通块（连通块就是前面说的“紧密排列”），的方案数。

新加入第 (i) 个数时，分三种情况：

可以让它自己作为一个连通块。此时它不会对 (j) 产生任何贡献，因为之后它两边的数一定都比它大。连通块数量 (k) 会加 (1)。
可以让它贴在某个连通块的一侧。此时它对 (j) 的贡献为 (i)。连通块数量 (k) 不变。
可以让它连接起两个连通块（也就是左右各紧挨着一个连通块）。此时它对 (j) 的贡献为 (2i)。连通块数量 (k) 减少 (1)。

注意，我们加入一个数时，并没有确定它在排列里的真实位置，而是确定它和其他已加入的数的相对位置关系。

这个 DP 的时间复杂度为 (O(n^4))，因为第 2 维 (j) 大小为 (n(n + 1))，即 (O(n^2)) 的。通过使用滚动数组，可以将空间复杂度优化到 (O(n^3))。

这个 DP 的结果就是 (dp[n][s][1]) ((s in[0,n(n + 1)]))。答案就是 (sum_{s = 0}^{n(n + 1)}dp[n][s][1] imes{L-s-1+nchoose n})。

问题转化为求 ({L-s-1+nchoose n})。因为 (m) 不一定是质数，我们不好求逆元。所以要用更奇妙的方法。

考虑第一维 (L-s-1+n) 的值，会有一个上下界，即 (l=max(L-n(n+1)-1+n,n),r=L-1+n)。两者相差是 (O(n^2)) 的。

先用矩阵快速幂求出 ({lchoose 0dots n})。再通过 ({ichoose j}={i-1choose j-1}+{i-1choose j}) 这个式子递推出所有 ({ldots rchoose 0dots n})。第一维大小 (O(n^2))，第二维大小 (O(n))。

时间复杂度 (O(n^3log L+n^3))。

总时间复杂度 (O(n^4+n^3log L))。

参考代码

// problem: LOJ6069
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fi first
#define se second
#define SZ(x) ((int)(x).size())

typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;

template<typename T> inline void ckmax(T& x, T y) { x = (y > x ? y : x); }
template<typename T> inline void ckmin(T& x, T y) { x = (y < x ? y : x); }

const int MAXN = 100;
int n, L, MOD, S;
int dp[2][MAXN * (MAXN + 1) + 5][MAXN + 5];

inline int mod1(int x) { return x < MOD ? x : x - MOD; }
inline int mod2(int x) { return x < 0 ? x + MOD : x; }
inline void add(int &x, int y) { x = mod1(x + y); }
inline void sub(int &x, int y) { x = mod2(x - y); }

struct Matrix {
	int a[MAXN + 5][MAXN + 5];
	void identity() {
		for(int i = 0; i <= n; ++i) {
			for(int j = 0; j <= n; ++j) {
				a[i][j] = (i == j);
			}
		}
	}
	Matrix() {
		memset(a, 0, sizeof(a));
	}
};
Matrix operator * (const Matrix& X, const Matrix& Y) {
	Matrix Z;
	for(int i = 0; i <= n; ++i) {
		for(int j = 0; j <= n; ++j) {
			for(int k = 0; k <= n; ++k) {
				Z.a[i][j] = ((ll)Z.a[i][j] + (ll)X.a[i][k] * Y.a[k][j]) % MOD;
			}
		}
	}
	return Z;
}
Matrix mat_pow(Matrix X, int i) {
	Matrix Y; Y.identity();
	while(i) {
		if(i & 1) Y = Y * X;
		X = X * X;
		i >>= 1;
	}
	return Y;
}

int comb[MAXN * (MAXN + 1) + 5][MAXN + 5];
void comb_init(int l, int r) {
	Matrix trans;
	trans.a[0][0] = 1;
	for(int i = 1; i <= n; ++i) {
		trans.a[i - 1][i] = 1;
		trans.a[i][i] = 1;
	}
	Matrix res = mat_pow(trans, l);
	for(int i = 0; i <= n; ++i) {
		comb[0][i] = res.a[0][i];
	}
	for(int i = l + 1; i <= r; ++i) {
		comb[i - l][0] = 1;
		for(int j = 1; j <= n; ++j) {
			comb[i - l][j] = mod1(comb[i - 1 - l][j - 1] + comb[i - 1 - l][j]);
		}
	}
}
int main() {
	cin >> n >> L >> MOD;
	dp[0][0][0] = 1;
	S = 0;
	for(int i = 1; i <= n; ++i) {
		int cur = (i & 1);
		int lst = (cur ^ 1);
		for(int j = 0; j <= S; ++j) for(int k = 0; k < i; ++k) dp[cur][j][k] = 0;
		for(int j = 0; j <= S; ++j) {
			for(int k = 0; k < i; ++k) if(dp[lst][j][k]) {
				add(dp[cur][j][k + 1], (ll)dp[lst][j][k] * (k + 1) % MOD);
				add(dp[cur][j + i][k], (ll)dp[lst][j][k] * 2 * k % MOD);
				if(k >= 2) add(dp[cur][j + i * 2][k - 1], (ll)dp[lst][j][k] * (k - 1) % MOD);
			}
		}
		S += i * 2;
	}
	// for(int i = 0; i <= S; ++i) cout << dp[n & 1][i][1] << " "; cout << endl;
	int l = max(L - S - 1 + n, n), r = L - 1 + n;
	if(l > r) { cout << 0 << endl; return 0; }
	comb_init(l, r);
	int ans = 0;
	for(int s = 0; s <= S; ++s) if(dp[n & 1][s][1]) {
		if(s + 1 > L) break;
		assert(L - s - 1 + n >= l && L - s - 1 + n <= r);
		ans = ((ll)ans + (ll)comb[L - s - 1 + n - l][n] * dp[n & 1][s][1]) % MOD;
	}
	cout << ans << endl;
	return 0;
}