行列
(math/matrix.hpp)

View this file on GitHub
Last update: 2022-10-20 00:25:59+09:00
Include: #include "math/matrix.hpp"

概要

集合 $V$ について定めた代数的構造上の $n \times m$ 行列に関する各種演算をサポートする。

詳細

<typename S> struct matrix
行列本体。
- typename S
  テンプレート引数として与える構造体。
  - using val_t
    $V$ を表す型。各種演算、比較、複合代入をサポートする必要がある。
  - val_t zero()
    加法単位元 $0$ を返す。
  - val_t one()
    乗法単位元 $1$ を返す。
- (constructor)(int n, int m)
  大きさ $n \times m$ 、全要素 $\mathrm{0}$ で初期化。
- (constructor)(vector<vector<V>> src)
  $src$ で初期化。
- 要素アクセス
  二次元配列のように、[]演算子を用いて行う。
- int height()
  高さ $n$ を返す。
- int width()
  幅 $m$ を返す。
- static matrix id(int n)
  $n \times n$ の単位行列を返す。
- S::val_t det()
  行列式を返す。定義に従った計算ではなく、上三角化を経由する関係上、S::val_t の除法を要請する。 $O(n^3)$ 時間。
- matrix inv()
  逆行列を返す。正則でない場合の動作は未定義。 $O(n^3)$ 時間。
- matrix pow(ll k)
  $k$ 乗した結果を返す。S::val_t が半環 (加法について可換モノイド、乗法についてモノイド、分配的、加法単位元は乗法吸収元) であることを要請する。 $O(n^3 \log k)$ 時間。
- 四則演算、比較、複合代入
struct double_field
実数体を載せるときにmatrixにテンプレート引数 S として与える。誤差周りの事情から、特殊化を行っている。

参考

うしさんのライブラリ

Depends on

template.hpp

Verified with

test/judge.yosupo.jp/Determinant_of_Matrix.0.test.cpp

Code

#pragma once

#include "../template.hpp"

template <typename S> struct matrix {
    using V = typename S::V;
    vector<vector<V>> val;
    matrix(int n, int m) : matrix(vector(n, vector(m, S::zero()))) {}
    matrix(const vector<vector<V>> &src) : val(src) {}
    vector<V> &operator[](int i) { return val[i]; }
    const vector<V> &operator[](int i) const { return val[i]; }
    int height() const { return val.size(); }
    int width() const { return val[0].size(); }
    static matrix id(int n) {
        matrix ret(n, n);
        for (int i : rep(n)) ret[i][i] = S::one();
        return ret;
    }
    void row_add(int i, int j, V a) {
        for (int k : rep(width())) { val[i][k] += val[j][k] * a; }
    }
    bool place_nonzero(int i, int j) {
        for (int k : rep(i, height())) {
            if (val[k][j] != S::zero()) {
                if (k > i) row_add(i, k, S::one());
                break;
            }
        }
        return val[i][j] != S::zero();
    }
    matrix upper_triangular() const {
        matrix ret(*this);
        for (int i = 0, j = 0; i < height() && j < width(); j++) {
            if (!ret.place_nonzero(i, j)) continue;
            for (int k : rep(i + 1, height())) ret.row_add(k, i, -ret[k][j] / ret[i][j]);
            i++;
        }
        return ret;
    }
    V det() const {
        V ret = S::one();
        matrix ut = upper_triangular();
        for (int i : rep(height())) ret *= ut[i][i];
        return ret;
    }
    matrix inv() const {
        matrix ex(height(), width() << 1);
        for (int i : rep(height())) {
            for (int j : rep(width())) { ex[i][j] = val[i][j]; }
        }
        for (int i : rep(height())) ex[i][width() + i] = S::one();
        matrix ut = ex.upper_triangular();
        for (int i : per(height())) {
            ut.row_add(i, i, S::one() / ut[i][i] - S::one());
            for (int j : rep(i)) ut.row_add(j, i, -ut[j][i] / ut[i][i]);
        }
        matrix ret(height(), width());
        for (int i : rep(height())) {
            for (int j : rep(width())) { ret[i][j] = ut[i][width() + j]; }
        }
        return ret;
    }
    matrix pow(ll k) const {
        matrix ret = matrix::id(height()), mul(*this);
        while (k) {
            if (k & 1) ret *= mul;
            mul *= mul;
            k >>= 1;
        }
        return ret;
    }
    matrix &operator+=(const matrix &a) {
        for (int i : rep(height())) {
            for (int j : rep(width())) { val[i][j] += a[i][j]; }
        }
        return *this;
    }
    matrix &operator-=(const matrix &a) {
        for (int i : rep(height())) {
            for (int j : rep(width())) { val[i][j] -= a[i][j]; }
        }
        return *this;
    }
    matrix &operator*=(const matrix &a) {
        matrix res(height(), a.width());
        for (int i : rep(height())) {
            for (int j : rep(a.width())) {
                for (int k : rep(width())) { res[i][j] += val[i][k] * a[k][j]; }
            }
        }
        val.swap(res.val);
        return *this;
    }
    matrix &operator/=(const matrix &a) { return *this *= a.inv(); }
    bool operator==(const matrix &a) const { return val == a.val; }
    bool operator!=(const matrix &a) const { return rel_ops::operator!=(*this, a); }
    matrix operator+() const { return *this; }
    matrix operator-() const { return matrix(height(), width()) -= *this; }
    matrix operator+(const matrix &a) const { return matrix(*this) += a; }
    matrix operator-(const matrix &a) const { return matrix(*this) -= a; }
    matrix operator*(const matrix &a) const { return matrix(*this) *= a; }
    matrix operator/(const matrix &a) const { return matrix(*this) /= a; }
};

struct double_field {
    using V = double;
    static V zero() { return 0.0; }
    static V one() { return 1.0; }
};

template <> bool matrix<double_field>::place_nonzero(int i, int j) {
    static constexpr double EPS = 1e-12;
    for (int k : rep(i + 1, height())) {
        if (abs(val[k][j]) > abs(val[i][j])) {
            swap(val[i], val[k]);
            row_add(i, i, -2.0);
        }
    }
    return abs(val[i][j]) > EPS;
}

#line 2 "math/matrix.hpp"

#line 2 "template.hpp"

#include <bits/stdc++.h>
using namespace std;

#define all(a) begin(a), end(a)
#define rall(a) rbegin(a), rend(a)
#define uniq(a) (a).erase(unique(all(a)), (a).end())
#define t0 first
#define t1 second
using ll = long long;
using ull = unsigned long long;
using pll = pair<ll, ll>;
using vll = vector<ll>;
constexpr double pi = 3.14159265358979323846;
constexpr ll dy[9] = {0, 1, 0, -1, 1, 1, -1, -1, 0};
constexpr ll dx[9] = {1, 0, -1, 0, 1, -1, -1, 1, 0};
constexpr ll sign(ll a) { return (a > 0) - (a < 0); }
constexpr ll fdiv(ll a, ll b) { return a / b - ((a ^ b) < 0 && a % b); }
constexpr ll cdiv(ll a, ll b) { return -fdiv(-a, b); }
constexpr ll pw(ll n) { return 1ll << n; }
constexpr ll flg(ll n) { return 63 - __builtin_clzll(n); }
constexpr ll clg(ll n) { return n == 1 ? 0 : flg(n - 1) + 1; }
constexpr ll safemod(ll x, ll mod) { return (x % mod + mod) % mod; }
template <typename T> using priority_queue_rev = priority_queue<T, vector<T>, greater<T>>;
template <typename T> constexpr T sq(const T &a) { return a * a; }
template <typename T, typename U> constexpr bool chmax(T &a, const U &b) { return a < b ? a = b, true : false; }
template <typename T, typename U> constexpr bool chmin(T &a, const U &b) { return a > b ? a = b, true : false; }
template <typename T, typename U> ostream &operator<<(ostream &os, const pair<T, U> &a) {
    os << "(" << a.first << ", " << a.second << ")";
    return os;
}
template <typename T, typename U, typename V> ostream &operator<<(ostream &os, const tuple<T, U, V> &a) {
    os << "(" << get<0>(a) << ", " << get<1>(a) << ", " << get<2>(a) << ")";
    return os;
}
template <typename T> ostream &operator<<(ostream &os, const vector<T> &a) {
    os << "(";
    for (auto itr = a.begin(); itr != a.end(); ++itr) os << *itr << (next(itr) != a.end() ? ", " : "");
    os << ")";
    return os;
}
template <typename T> ostream &operator<<(ostream &os, const set<T> &a) {
    os << "(";
    for (auto itr = a.begin(); itr != a.end(); ++itr) os << *itr << (next(itr) != a.end() ? ", " : "");
    os << ")";
    return os;
}
template <typename T> ostream &operator<<(ostream &os, const multiset<T> &a) {
    os << "(";
    for (auto itr = a.begin(); itr != a.end(); ++itr) os << *itr << (next(itr) != a.end() ? ", " : "");
    os << ")";
    return os;
}
template <typename T, typename U> ostream &operator<<(ostream &os, const map<T, U> &a) {
    os << "(";
    for (auto itr = a.begin(); itr != a.end(); ++itr) os << *itr << (next(itr) != a.end() ? ", " : "");
    os << ")";
    return os;
}
#ifdef ONLINE_JUDGE
#define dump(...) (void(0))
#else
void debug() { cerr << endl; }
template <typename Head, typename... Tail> void debug(Head &&head, Tail &&... tail) {
    cerr << head;
    if (sizeof...(Tail)) cerr << ", ";
    debug(tail...);
}
#define dump(...) cerr << __LINE__ << ": " << #__VA_ARGS__ << " = ", debug(__VA_ARGS__)
#endif
struct rep {
    struct itr {
        ll v;
        itr(ll v) : v(v) {}
        void operator++() { ++v; }
        ll operator*() const { return v; }
        bool operator!=(itr i) const { return v < *i; }
    };
    ll l, r;
    rep(ll l, ll r) : l(l), r(r) {}
    rep(ll r) : rep(0, r) {}
    itr begin() const { return l; };
    itr end() const { return r; };
};
struct per {
    struct itr {
        ll v;
        itr(ll v) : v(v) {}
        void operator++() { --v; }
        ll operator*() const { return v; }
        bool operator!=(itr i) const { return v > *i; }
    };
    ll l, r;
    per(ll l, ll r) : l(l), r(r) {}
    per(ll r) : per(0, r) {}
    itr begin() const { return r - 1; };
    itr end() const { return l - 1; };
};
struct io_setup {
    static constexpr int PREC = 20;
    io_setup() {
        cout << fixed << setprecision(PREC);
        cerr << fixed << setprecision(PREC);
    };
} iOS;
#line 4 "math/matrix.hpp"

template <typename S> struct matrix {
    using V = typename S::V;
    vector<vector<V>> val;
    matrix(int n, int m) : matrix(vector(n, vector(m, S::zero()))) {}
    matrix(const vector<vector<V>> &src) : val(src) {}
    vector<V> &operator[](int i) { return val[i]; }
    const vector<V> &operator[](int i) const { return val[i]; }
    int height() const { return val.size(); }
    int width() const { return val[0].size(); }
    static matrix id(int n) {
        matrix ret(n, n);
        for (int i : rep(n)) ret[i][i] = S::one();
        return ret;
    }
    void row_add(int i, int j, V a) {
        for (int k : rep(width())) { val[i][k] += val[j][k] * a; }
    }
    bool place_nonzero(int i, int j) {
        for (int k : rep(i, height())) {
            if (val[k][j] != S::zero()) {
                if (k > i) row_add(i, k, S::one());
                break;
            }
        }
        return val[i][j] != S::zero();
    }
    matrix upper_triangular() const {
        matrix ret(*this);
        for (int i = 0, j = 0; i < height() && j < width(); j++) {
            if (!ret.place_nonzero(i, j)) continue;
            for (int k : rep(i + 1, height())) ret.row_add(k, i, -ret[k][j] / ret[i][j]);
            i++;
        }
        return ret;
    }
    V det() const {
        V ret = S::one();
        matrix ut = upper_triangular();
        for (int i : rep(height())) ret *= ut[i][i];
        return ret;
    }
    matrix inv() const {
        matrix ex(height(), width() << 1);
        for (int i : rep(height())) {
            for (int j : rep(width())) { ex[i][j] = val[i][j]; }
        }
        for (int i : rep(height())) ex[i][width() + i] = S::one();
        matrix ut = ex.upper_triangular();
        for (int i : per(height())) {
            ut.row_add(i, i, S::one() / ut[i][i] - S::one());
            for (int j : rep(i)) ut.row_add(j, i, -ut[j][i] / ut[i][i]);
        }
        matrix ret(height(), width());
        for (int i : rep(height())) {
            for (int j : rep(width())) { ret[i][j] = ut[i][width() + j]; }
        }
        return ret;
    }
    matrix pow(ll k) const {
        matrix ret = matrix::id(height()), mul(*this);
        while (k) {
            if (k & 1) ret *= mul;
            mul *= mul;
            k >>= 1;
        }
        return ret;
    }
    matrix &operator+=(const matrix &a) {
        for (int i : rep(height())) {
            for (int j : rep(width())) { val[i][j] += a[i][j]; }
        }
        return *this;
    }
    matrix &operator-=(const matrix &a) {
        for (int i : rep(height())) {
            for (int j : rep(width())) { val[i][j] -= a[i][j]; }
        }
        return *this;
    }
    matrix &operator*=(const matrix &a) {
        matrix res(height(), a.width());
        for (int i : rep(height())) {
            for (int j : rep(a.width())) {
                for (int k : rep(width())) { res[i][j] += val[i][k] * a[k][j]; }
            }
        }
        val.swap(res.val);
        return *this;
    }
    matrix &operator/=(const matrix &a) { return *this *= a.inv(); }
    bool operator==(const matrix &a) const { return val == a.val; }
    bool operator!=(const matrix &a) const { return rel_ops::operator!=(*this, a); }
    matrix operator+() const { return *this; }
    matrix operator-() const { return matrix(height(), width()) -= *this; }
    matrix operator+(const matrix &a) const { return matrix(*this) += a; }
    matrix operator-(const matrix &a) const { return matrix(*this) -= a; }
    matrix operator*(const matrix &a) const { return matrix(*this) *= a; }
    matrix operator/(const matrix &a) const { return matrix(*this) /= a; }
};

struct double_field {
    using V = double;
    static V zero() { return 0.0; }
    static V one() { return 1.0; }
};

template <> bool matrix<double_field>::place_nonzero(int i, int j) {
    static constexpr double EPS = 1e-12;
    for (int k : rep(i + 1, height())) {
        if (abs(val[k][j]) > abs(val[i][j])) {
            swap(val[i], val[k]);
            row_add(i, i, -2.0);
        }
    }
    return abs(val[i][j]) > EPS;
}

行列 (math/matrix.hpp)

概要

詳細

参考

Depends on

Verified with

Code

行列
(math/matrix.hpp)