P3390 【模板】矩阵快速幂 · znzryb
## 思路讲解
```cpp
// 矩阵快速幂：必须是方阵
template <class T>
Matrix<T> mat_pow(Matrix<T> base, long long exp) {
assert(base.n == base.m);
Matrix<T> res = Matrix<T>::identity(base.n);
while (exp > 0) {
if (exp & 1) res = res * base;
base = base * base;
exp >>= 1;
}
return res;
}
```
## AC代码
::: 源代码
```cpp
// by Gospel_rock
// title: P3390 【模板】矩阵快速幂
// GROUP_NAME: // ONLINE_JUDGE_NAME: Luogu
#include <bits/stdc++.h>
#define all(vec) vec.begin(),vec.end()
#define PB push_back
#define PF push_front
#define EB emplace_back
#define EF emplace_front
#define EM emplace
#define FI first
#define SE second
#define pct __builtin_popcountll
#define ctz __builtin_ctzll
#define SZ(a) ((long long) a.size())
#define FOR(i, a, b) for (int i = (a); i <= (b); ++i)
#define ROF(i, a, b) for (int i = (a); i >= (b); --i)
#define debug(var) cerr << #var <<":"<<var<<"\n";
#define cend cerr<<"\n-----------\n"
#define fsp(x) fixed<<setprecision(x)
using namespace std;
using ll = long long;using ull = unsigned long long;
using DB=double;using LD=long double;
// using i128 = __int128;
using CD=complex<double>;
static constexpr ll MAXN=(ll)1e6+10,INF=(1ll<<61)-1;
static constexpr ll mod=(ll)1e9+7; // 998244353;
static constexpr double eps=1e-8;const double pi=acos(-1.0);
ll lT;
// T: 底层存储类型 (uint32_t / uint64_t)// MOD: 模数
// Wide: 乘法时用的更宽类型
template <class T, T MOD, class Wide = ull>
struct modint {
T v;
modint(ll x = 0) {
ll t = x % (ll)MOD;
if (t < 0) t += MOD;
v = (T)t;
}
T val() const { return v; }
modint& operator+=(const modint &o) {
T x = v + o.v;
if (x >= MOD) x -= MOD;
v = x;
return *this;
}
modint& operator-=(const modint &o) {
T x = v >= o.v ? v - o.v : (T)(v + MOD - o.v);
v = x;
return *this;
}
modint& operator*=(const modint &o) {
Wide t = (Wide)v * (Wide)o.v;
v = (T)(t % MOD);
return *this;
}
// 幂
static modint powmod(modint a, ll e) {
modint r = 1;
while (e) {
if (e & 1) r *= a;
a *= a;
e >>= 1;
}
return r;
}
// 逆元（假定 MOD 是质数）
static modint inv(modint a) {
return powmod(a, (ll)MOD - 2);
}
modint& operator/=(const modint &o) {
return *this *= inv(o);
}
// 一元
modint operator+() const { return *this; }
modint operator-() const { return modint(0) - *this; }
// 二元
friend modint operator+(modint a, const modint &b) { return a += b; }
friend modint operator-(modint a, const modint &b) { return a -= b; }
friend modint operator*(modint a, const modint &b) { return a *= b; }
friend modint operator/(modint a, const modint &b) { return a /= b; }
friend bool operator==(const modint &a, const modint &b) { return a.v == b.v; }
friend bool operator!=(const modint &a, const modint &b) { return a.v != b.v; }
friend ostream& operator<<(ostream &os, const modint &x) {
return os << x.v;
}
friend istream& operator>>(istream &is, modint &x) {
ll t;
if (!(is >> t)) return is;
x = modint(t);
return is;
}
};
using modi=modint<uint32_t,mod>;
// 2025/11/5 AC https://www.luogu.com.cn/record/245401024
template <class T>
struct Matrix {
int n, m;
vector<vector<T>> a; // 1-based
Matrix(): n(0), m(0) {}
Matrix(int n_, int m_, const T &val = T()) {
assign(n_, m_, val);
}
// 新增：直接用 1-based 的 vv 构造
Matrix(const vector<vector<T>> &vv) {
load(vv);
}
// 你已经有一个 1-based 的 vv（vv.size() == n+1, vv[1].size() == m+1）
void load(const vector<vector<T>> &vv) {
n = (int)vv.size() - 1;
m = n ? (int)vv[1].size() - 1 : 0;
a = vv; // 直接拷
}
void assign(int n_, int m_, const T &val = T()) {
n = n_;
m = m_;
a.assign(n + 1, vector<T>(m + 1, val));
}
vector<T>& operator[](int i) {
return a[i]; // A[i][j] -> (A.operator[](i))[j]
}
// 在const的时候调用
const vector<T>& operator[](int i) const {
return a[i];
}
Matrix operator+(const Matrix &rhs) const {
assert(n == rhs.n && m == rhs.m);
Matrix res(n, m);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m; ++j)
res.a[i][j] = a[i][j] + rhs.a[i][j];
return res;
}
Matrix& operator+=(const Matrix &rhs) {
return *this = *this + rhs;
}
Matrix operator-(const Matrix &rhs) const {
assert(n == rhs.n && m == rhs.m);
Matrix res(n, m);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m; ++j)
res.a[i][j] = a[i][j] - rhs.a[i][j];
return res;
}
Matrix& operator-=(const Matrix &rhs) {
return *this = *this - rhs;
}
Matrix operator*(const Matrix &rhs) const {
assert(m == rhs.n);
Matrix res(n, rhs.m, T());
for (int i = 1; i <= n; ++i) {
for (int k = 1; k <= m; ++k) {
if (a[i][k] == T()) continue;
for (int j = 1; j <= rhs.m; ++j) {
res.a[i][j] = res.a[i][j] + a[i][k] * rhs.a[k][j];
}
}
}
return res;
}
Matrix& operator*=(const Matrix &rhs) {
return *this = *this * rhs;
}
static Matrix identity(int n) {
Matrix I(n, n, T());
FOR(i,1,n) I.a[i][i] = T(1);
return I;
}
friend ostream& operator<<(ostream &os,const Matrix& ma) {
FOR(i,1,ma.n) {
FOR(j,1,ma.m) {
os<<ma[i][j]<<" ";
}
os<<"\n";
}
return os;
}
};
// 矩阵快速幂：必须是方阵
template <class T>
Matrix<T> mat_pow(Matrix<T> base, long long exp) {
assert(base.n == base.m);
Matrix<T> res = Matrix<T>::identity(base.n);
while (exp > 0) {
if (exp & 1) res = res * base;
base = base * base;
exp >>= 1;
}
return res;
}
void Solve(){
ll n,k;
cin>>n>>k;
vector<vector<modi> > A(n+1,vector<modi>(n+1));
FOR(i,1,n) {
FOR(j,1,n) {
cin>>A[i][j];
}
}
Matrix mat(A);
mat=mat_pow(mat,k);
cout<<mat;
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
// cin>>lT;
// while(lT--) Solve();
return 0;
}
/*
*/
```
:::
## 心路历程（WA，TLE，MLE……）