题目来源：LibreOJ #6069. 「2017 山东一轮集训 Day4」塔
题目大意
我们的主人公 yzh 有 \(L\) 个妹子，她们站成一排，依次编号为 \(1\dots 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\) 取模。
数据范围 \(1\leq L\leq 10^9\)，\(1\leq n\leq 100\)，\(1\leq m\leq 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 + n\choose n}\)。
\(s\) 最大不超过 \(\sum_{i=1}^{n}2i\)，即 \(n(n+1)\)。
考虑先对每个 \(s\in [0,n(n+1)]\)，求出 \(\sum_{i = 2}^{n}\max(p_i,p_{i - 1}) = s\) 的排列 \(p\) 的数量。
用 DP 求。设 \(dp[i][j][k]\) 表示从小到大考虑了 \(1\dots 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]\times{L-s-1+n\choose n}\)。
问题转化为求 \({L-s-1+n\choose n}\)。因为 \(m\) 不一定是质数，我们不好求逆元。所以要用更奇妙的方法。
考虑第一维 \(L-s-1+n\) 的值，会有一个上下界，即 \(l=\max(L-n(n+1)-1+n,n),r=L-1+n\)。
先用矩阵快速幂求出 \({l\choose 0\dots n}\)。再通过 \({i\choose j}={i-1\choose j-1}+{i-1\choose j}\) 这个式子递推出所有 \({l\dots r\choose 0\dots n}\)。
时间复杂度 \(O(n^3\log L+n^3)\)。
总时间复杂度 \(O(n^4+n^3\log 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;
}