P15569 [COCI 2025/2026 #5] 结构 / Struktura 题解
tag : 矩阵快速幂 矩阵加速 DP 转移
赛时这道题写了 \(O(n)\) DP,得了 76 分,但不太熟悉矩阵加速,没有写出正解
这道题有更简单的方法,但为了练习矩阵加速 DP,用了如下方法
题意
给定 \(n,k\),独立、等概率地从 \(1\) 到 \(k\) 中取 \(n\) 个数,得到数组 \(a_1,\dots,a_n\)。称 \(a\) 是结构,当且仅当同时满足
- \(1,2,\dots,n\) 在 \(a\) 中各出现恰好一次,即 \(a\) 是 \(1\) 到 \(n\) 的一个排列;
- 对每个 \(i\),\(|a_i+i-n-1|\le 1\)。
求随机数组恰为结构的概率,模 \(10^9+7\) 输出,即
数据范围 \(n,k\le 10^9\)。
第二条等价于每个位置只有三种候选值
落出 \([1,n]\) 的候选值不合法。另外 \(k 时答案为 \(0\),因为排列里的 \(n\) 也必须取自 \(1\) 到 \(k\)。
显然,这道题有一个比较简单的暴力 DP 转移
显然,每个 \(a_i\) 的取值只有三种,从 \(i=1\) 到 \(n\) 分别如下:
| \(i\) | \(1\) | \(2\) | \(3\) | \(\cdots\) | \(n-2\) | \(n-1\) | \(n\) |
|---|---|---|---|---|---|---|---|
| ① | \(n+1\) | \(n\) | \(n-1\) | \(\cdots\) | \(4\) | \(3\) | \(2\) |
| ② | \(n\) | \(n-1\) | \(n-2\) | \(\cdots\) | \(3\) | \(2\) | \(1\) |
| ③ | \(n-1\) | \(n-2\) | \(n-3\) | \(\cdots\) | \(2\) | \(1\) | \(0\) |
\(f[i][a][b]\) : 表示第 \(i\) 位选第 \(a\) 个,第 \(i - 1\) 位选第 \(b\) 个,注意边界 \(1\) 到 \(n\)
// init
f[2][1][1] = 0, f[2][1][3] = 1;
f[2][2][1] = 0, f[2][2][2] = 1;
f[2][3][1] = 0, f[2][3][2] = 1, f[2][3][3] = 1;for (int i = 2; i < n; i ++) {f[i + 1][1][1] = (f[i + 1][1][1] + f[i][1][1]) % MOD;f[i + 1][1][3] = (f[i + 1][1][3] + f[i][3][1] + f[i][3][2]) % MOD;f[i + 1][2][1] = (f[i + 1][2][1] + f[i][1][1] + f[i][1][3]) % MOD;f[i + 1][2][2] = (f[i + 1][2][2] + f[i][2][1] + f[i][2][2]) % MOD;f[i + 1][3][1] = (f[i + 1][3][1] + f[i][1][1] + f[i][1][3]) % MOD;f[i + 1][3][2] = (f[i + 1][3][2] + f[i][2][1] + f[i][2][2]) % MOD;f[i + 1][3][3] = (f[i + 1][3][3] + f[i][3][1] + f[i][3][2] + f[i][3][3]) % MOD;
}cnt = (cnt + f[n][1][1] + f[n][1][3]) % MOD;
cnt = (cnt + f[n][2][1] + f[n][2][2]) % MOD;
能够用矩阵加速的五个条件(\(AI\)总结)
- 一步能写成加权和 :新状态里每一项都是旧状态的“常数 × 旧值”再相加,比如 \(f i = f i − 1 + f i − 2 、 b 3 = a 7 + a 8\) 。不能出现两个状态相乘、开平方、绝对值、条件判断。
- 权重不随步数变 : \(f i = f i − 1 + f i − 2\) 里系数是 1 , 1 ,每步都一样;如果是 \(f i = i ⋅ f i − 1\) 或 \(f i = 2 i f i − 1\) ,系数每步都变,就不能直接用。
- 状态个数固定 :每一步只需要固定几个数(比如本题 9 个、斐波那契 2 个)。如果需要记住整段历史、或者状态数是 \(2^n\) 这种规模,就不行。
- 只要“某一时刻的状态” :答案能直接从第 \(n\) 步的状态里取出来。
- 再加一个现实条件: \(n\) 要足够大 ,否则直接循环反而更快。
所以,显然本题可以用矩阵加速
矩阵加速DP转移:
令 \(id ( a , b ) = ( a − 1 ) ⋅ 3 + b\)
所以可以把暴力 DP 转移写成如下形式 (重要\(! ! !\))
其中
考虑如何推出 \(base\) 矩阵
由暴力 DP
\(f_{i,1}=f_{i-1,1}\)
\(f_{i,3}=f_{i-1,7}+f_{i-1,8}\)
\(f_{i,4}=f_{i-1,1}+f_{i-1,3}\)
\(f_{i,5}=f_{i-1,4}+f_{i-1,5}\)
\(f_{i,7}=f_{i-1,1}+f_{i-1,3}\)
\(f_{i,8}=f_{i-1,4}+f_{i-1,5}\)
\(f_{i,9}=f_{i-1,7}+f_{i-1,8}+f_{i-1,9}\)
对于 \(f_{i,1}\) 由\(f_{i,1}=f_{i-1,1}\) 可以得到:
同理可以得到整个 \(base\) 矩阵
于是 \(b[j] = \sum_{i=1}^{9} a[i] * base[i][j]\)
于是直接开两个数组(无需用矩阵存\(b\), \(a\))
Matrix p = qpow(base, n - 2);a[3] = a[5] = a[8] = a[9] = 1;
for (int j = 1; j <= 9; j ++) {for (int i = 1; i <= 9; i ++) {b[j] = (b[j] + a[i] * p.mat[i][j]) % MOD;}
}
答案:b[1] + b[3] + b[4] + b[5]
AC code
#include
using namespace std;using ll = long long;const int N = 10;
const ll MOD = 1e9 + 7;ll nn, k;ll a[N], b[N];struct Matrix {int n;ll mat[N][N];Matrix(int _n = 0) {n = _n;memset(mat, 0, sizeof(mat));}void E() {for (int i = 1; i <= n; i ++) mat[i][i] = 1;}
} ;Matrix mul(Matrix a, Matrix b) {int n = a.n;Matrix ret(n);for (int i = 1; i <= n; i ++) {for (int j = 1; j <= n; j ++) {for (int k = 1; k <= n; k ++) {ret.mat[i][j] = (ret.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % MOD;}}}return ret;
}Matrix qpow(Matrix a, ll b) {int n = a.n;Matrix res(n);res.E();while (b) {if (b & 1) res = mul(res, a);a = mul(a, a);b >>= 1;}return res;
}ll q_pow(ll a, ll b) {ll res = 1;while (b) {if (b & 1) res = res * a % MOD;a = a * a % MOD;b >>= 1;}return res;
}int main () {freopen("struktura.in", "r", stdin);freopen("struktura.out", "w", stdout);ios::sync_with_stdio(false);cin.tie(0);cin >> nn >> k;if (nn == 1) {cout << q_pow(k, MOD - 2) << '\n';return 0;}if (k < nn) {cout << 0 << '\n';return 0;}Matrix base(9);base.mat[1][1] = base.mat[1][4] = base.mat[1][7] =base.mat[3][4] = base.mat[3][7] =base.mat[4][5] = base.mat[4][8] = base.mat[5][5] = base.mat[5][8] = base.mat[7][3] = base.mat[7][9] = base.mat[8][3] = base.mat[8][9] = base.mat[9][9] = 1;Matrix p = qpow(base, nn - 2);a[3] = a[5] = a[8] = a[9] = 1;for (int j = 1; j <= 9; j ++) {for (int i = 1; i <= 9; i ++) {b[j] = (b[j] + a[i] * p.mat[i][j]) % MOD;}}ll cnt = (b[1] + b[3] + b[4] + b[5]) % MOD;ll tot = q_pow(k, nn) % MOD;cout << cnt * q_pow(tot, MOD - 2) % MOD << '\n';return 0;
}
注意特判 \(n == 1\)
务必注意矩阵的 \(i, j\),非常容易搞反