|
| 1 | +// 矩阵乘法 C = A * B |
| 2 | +vector<vector<ll>> matMul(const vector<vector<ll>>& A, const vector<vector<ll>>& B) { |
| 3 | + int n = A.size(); |
| 4 | + int m = B[0].size(); |
| 5 | + int K = B.size(); // A: n x K, B: K x m |
| 6 | + vector<vector<ll>> C(n, vector<ll>(m, 0)); |
| 7 | + for (int i = 0; i < n; ++i) { |
| 8 | + for (int k = 0; k < K; ++k) { |
| 9 | + ll aik = A[i][k]; |
| 10 | + if (aik == 0) continue; |
| 11 | + for (int j = 0; j < m; ++j) { |
| 12 | + C[i][j] += (aik * B[k][j]) % MOD; |
| 13 | + if (C[i][j] >= MOD) C[i][j] -= MOD; |
| 14 | + } |
| 15 | + } |
| 16 | + } |
| 17 | + return C; |
| 18 | +} |
| 19 | + |
| 20 | +// 矩阵自乘 (square) |
| 21 | +vector<vector<ll>> matMulSquare(const vector<vector<ll>>& A, const vector<vector<ll>>& B) { |
| 22 | + return matMul(A, B); |
| 23 | +} |
| 24 | + |
| 25 | +// 矩阵快速幂 T^e |
| 26 | +vector<vector<ll>> matPow(vector<vector<ll>> base, long long e) { |
| 27 | + int K = base.size(); |
| 28 | + // initialize res = identity |
| 29 | + vector<vector<ll>> res(K, vector<ll>(K, 0)); |
| 30 | + for (int i = 0; i < K; ++i) res[i][i] = 1; |
| 31 | + while (e > 0) { |
| 32 | + if (e & 1) res = matMul(res, base); |
| 33 | + base = matMul(base, base); |
| 34 | + e >>= 1; |
| 35 | + } |
| 36 | + return res; |
| 37 | +} |
| 38 | + |
| 39 | +// 矩阵乘向量 v' = M * v |
| 40 | +vector<ll> matVecMul(const vector<vector<ll>>& M, const vector<ll>& v) { |
| 41 | + int n = M.size(); |
| 42 | + int m = v.size(); |
| 43 | + vector<ll> r(n, 0); |
| 44 | + for (int i = 0; i < n; ++i) { |
| 45 | + ll sum = 0; |
| 46 | + for (int j = 0; j < m; ++j) { |
| 47 | + if (M[i][j] == 0) continue; |
| 48 | + sum += (M[i][j] * v[j]) % MOD; |
| 49 | + if (sum >= MOD) sum -= MOD; |
| 50 | + } |
| 51 | + r[i] = sum % MOD; |
| 52 | + } |
| 53 | + return r; |
| 54 | +} |
| 55 | + |
| 56 | +v_0 = {....}; |
| 57 | +T = {{...}}; |
| 58 | +vector<vector<ll>>T_p = matPow(T, n-1); |
| 59 | +vector<ll>v_n = matVecMul(T_p, v_0); |
0 commit comments