example/matrix.example.cpp

View this file on GitHub
Last update: 2025-01-10 23:29:41+09:00

Depends on

Code

#include "../Utility/template.hpp"
#include "../Utility/modint.hpp"
#include "../Math/matrix.hpp"
using mint = modint998244353;

int main() {
    int n = 3; 
    Matrix<mint>  mat(n, n, 3); // n * nで初期値が mint(0) の行列を生成

    mat[0][0] = 1; // (0, 0)成分を 1 に変更
    mat[0][1] += 10; //(0, 1)成分に 10　加算
    //mat[3][3] = 100; //配列外参照 assert 無し

    auto pow_mat = mat.pow(5); // matの 100 乗を受け取る。 matは破壊されない。

    auto print_vec = [](vec<vec<mint>> v) {
        rep(i, 0, v.size()) {
            rep(j, 0, v[i].size()) cout << v[i][j] << " ";
            cout << endl;
        }
    };

    print_vec(mat.dump());// dumpで返された二次元配列を出力
    print_vec(pow_mat.dump());// 5乗されている。

    Matrix<mint> v(n, 1, 1); //n * 1の行列。擬似的なベクトルとして扱える。

    Matrix<mint> res = mat * v;//行列同士の掛け算。行列累乗の際、このような処理を書くかもしれない。
    print_vec(res.dump());
    cout << res[0][0] << endl;

}

#line 1 "example/matrix.example.cpp"

#line 1 "Utility/template.hpp"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
#define rep(i, s, t) for (ll i = s; i < (ll)(t); i++)
#define rrep(i, s, t) for (ll i = (ll)(t) - 1; i >= (ll)(s); i--)
#define all(x) begin(x), end(x)

#define TT template <typename T>
TT using vec = vector<T>;
template <class T1, class T2> bool chmin(T1 &x, T2 y) {
    return x > y ? (x = y, true) : false;
}
template <class T1, class T2> bool chmax(T1 &x, T2 y) {
    return x < y ? (x = y, true) : false;
}
struct io_setup {
    io_setup() {
        ios::sync_with_stdio(false);
        std::cin.tie(nullptr);
        cout << fixed << setprecision(15);
    }
} io_setup;

/*
@brief verify用テンプレート
*/
#line 1 "Utility/modint.hpp"

// 動的mod : template<int mod> を消して、上の方で変数modを宣言
template <uint32_t mod> struct modint {
    using mm = modint;
    uint32_t x;
    modint() : x(0) {}
    TT modint(T a = 0) : x((ll(a) % mod + mod)) {
        if (x >= mod) x -= mod;
    }

    friend mm operator+(mm a, mm b) {
        a.x += b.x;
        if (a.x >= mod) a.x -= mod;
        return a;
    }
    friend mm operator-(mm a, mm b) {
        a.x -= b.x;
        if (a.x >= mod) a.x += mod;
        return a;
    }

    mm operator-() const { return mod - x; }

    //+と-だけで十分な場合、以下は省略して良いです。

    friend mm operator*(mm a, mm b) { return (uint64_t)(a.x) * b.x; }
    friend mm operator/(mm a, mm b) { return a * b.inv(); }
    friend mm &operator+=(mm &a, mm b) { return a = a + b; }
    friend mm &operator-=(mm &a, mm b) { return a = a - b; }
    friend mm &operator*=(mm &a, mm b) { return a = a * b; }
    friend mm &operator/=(mm &a, mm b) { return a = a * b.inv(); }

    mm inv() const {
        assert(x != 0);
        return pow(mod - 2);
    }
    mm pow(ll y) const {
        mm res = 1;
        mm v = *this;
        while (y) {
            if (y & 1) res *= v;
            v *= v;
            y /= 2;
        }
        return res;
    }

    friend istream &operator>>(istream &is, mm &a) {
        ll t;
        cin >> t;
        a = mm(t);
        return is;
    }

    friend ostream &operator<<(ostream &os, mm a) { return os << a.x; }

    bool operator==(mm a) { return x == a.x; }
    bool operator!=(mm a) { return x != a.x; }

    bool operator<(const mm &a) const { return x < a.x; }
};
using modint998244353 = modint<998244353>;
using modint1000000007 = modint<1'000'000'007>;
/*
@brief modint
*/
#line 1 "Math/matrix.hpp"
template <typename T> struct Matrix {
    int h, w;
    vector<vector<T>> d;
    Matrix() {}
    Matrix(int h, int w, T val = 0) : h(h), w(w), d(h, vector<T>(w, val)) {}
    Matrix(vector<vector<T>> const &dat) : h(dat.size()), w(0), d(dat) {
        if (h > 0) w = d[0].size();
    }

    static Matrix unit(int n) {
        Matrix uni(n, n, 0);
        for (int i = 0; i < n; i++) {
            uni[i][i] = 1;
        }
        return uni;
    }
    const vector<T> &operator[](int i) const { return d[i]; }
    vector<T> &operator[](int i) { return d[i]; }
    Matrix &operator*=(const Matrix &a) { return *this = (*this) * a; }
    Matrix operator*(const Matrix &a) const {
        assert(w == a.h);
        Matrix r(h, a.w);
        for (int i = 0; i < h; i++)
            for (int k = 0; k < w; k++)
                for (int j = 0; j < a.w; j++) {
                    r[i][j] += d[i][k] * a[k][j];
                }

        return r;
    }

    Matrix pow(ll t) const {
        assert(h == w);
        Matrix res = Matrix::unit(h);
        Matrix x = (*this);
        while (t > 0) {
            if (t & 1) res = res * x;
            x = x * x;
            t >>= 1;
        }
        return res;
    }

    tuple<Matrix, T, ll> gaussian_elimination(int w_limit = -1) const {
        if (w_limit == -1) w_limit = w;
        T k = 1;
        Matrix A = *this;
        int i1 = 0;
        for (int j = 0; j < w_limit; j++) {
            if (i1 >= h) break;
            for (int i2 = i1; i2 < h; i2++) {
                if (A[i2][j] != 0) {
                    swap(A[i1], A[i2]);
                    if (i1 != i2) k = -k;
                    break;
                }
            }
            if (A[i1][j] == 0) {
                continue;
            }
            T inv = 1 / A[i1][j];
            k *= A[i1][j];
            for (int jj = 0; jj < w; jj++) {
                A[i1][jj] *= inv;
            }
            for (int i = 0; i < h; i++)
                if (A[i][j] != 0 && i != i1) {
                    T c = -A[i][j];
                    for (int jj = 0; jj < w; jj++) {
                        A[i][jj] += A[i1][jj] * c;
                    }
                }
            i1++;
        }
        return make_tuple(A, k, i1);
    }

    ll rank() const {
        auto [dat, k, rnk] = (*this).gaussian_elimination();
        return rnk;
    }

    pair<vector<T>, bool> linear_equations() const {
        assert(h == w - 1);
        vector<T> ret(w - 1);
        auto [dat, p, rnk] = (*this).gaussian_elimination(w - 1);
        if (rnk != w - 1) return make_pair(ret, false);
        for (int i = 0; i < h; i++) {
            ret[i] = dat[i][w - 1];
        }
        return make_pair(ret, true);
    }

    pair<Matrix, bool> inv() const {
        assert(h == w);
        Matrix slv(h, w * 2);
        for (int i = 0; i < h; i++)
            for (int j = 0; j < w; j++) {
                slv[i][j] = (*this)[i][j];
            }
        for (int i = 0; i < h; i++) {
            slv[i][i + w] = 1;
        }

        auto [dat, p, rnk] = slv.gaussian_elimination(w);
        auto ret = Matrix::unit(h);
        if (rnk != h) return make_pair(ret, false);
        for (int i = 0; i < h; i++) {
            for (int j = 0; j < w; j++) {
                ret[i][j] = dat[i][j + w];
            }
        }
        return make_pair(ret, true);
    }

    T det() const {
        assert(h == w);
        auto [A, p, rnk] = (*this).gaussian_elimination();
        for (int i = 0; i < h; i++) p *= A[i][i];
        return p;
    }

    friend ostream &operator<<(ostream &os, Matrix a) {
        for (int i = 0; i < a.h; i++) {
            for (int j = 0; j < a.w; j++) {
                os << a[i][j] << (j != a.w - 1 ? " " : "");
            }
            os << (i != a.h - 1 ? "\n" : "");
        }
        return os;
    }
};
#line 5 "example/matrix.example.cpp"
using mint = modint998244353;

int main() {
    int n = 3; 
    Matrix<mint>  mat(n, n, 3); // n * nで初期値が mint(0) の行列を生成

    mat[0][0] = 1; // (0, 0)成分を 1 に変更
    mat[0][1] += 10; //(0, 1)成分に 10　加算
    //mat[3][3] = 100; //配列外参照 assert 無し

    auto pow_mat = mat.pow(5); // matの 100 乗を受け取る。 matは破壊されない。

    auto print_vec = [](vec<vec<mint>> v) {
        rep(i, 0, v.size()) {
            rep(j, 0, v[i].size()) cout << v[i][j] << " ";
            cout << endl;
        }
    };

    print_vec(mat.dump());// dumpで返された二次元配列を出力
    print_vec(pow_mat.dump());// 5乗されている。

    Matrix<mint> v(n, 1, 1); //n * 1の行列。擬似的なベクトルとして扱える。

    Matrix<mint> res = mat * v;//行列同士の掛け算。行列累乗の際、このような処理を書くかもしれない。
    print_vec(res.dump());
    cout << res[0][0] << endl;

}