題目來源：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\) 的數量。注意，\(n(n + 1)\) 只是我們粗略估計的 \(s\) 的一個上界，事實是達不到這個上界的，不過不要緊，不存在這樣的 \(s\) 時算方案數為 \(0\) 即可。
用 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\)。兩者相差是 \(O(n^2)\) 的。
先用矩陣快速冪求出 \({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^2)\)，第二維大小 \(O(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;
}