P3390 【模板】矩阵快速幂

思路讲解

// 矩阵快速幂：必须是方阵
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代码

源代码

// 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……）