結果

問題 No.2441 行列累乗
ユーザー Gandalfr
提出日時 2023-09-02 08:38:53
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 3 ms / 2,000 ms
コード長 10,398 bytes
コンパイル時間 2,166 ms
コンパイル使用メモリ 198,744 KB
最終ジャッジ日時 2025-02-16 17:56:13
ジャッジサーバーID
（参考情報） 		judge4 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 20
権限があれば一括ダウンロードができます

ソースコード

diff #
プレゼンテーションモードにする

#line 1 "playspace/main.cpp"
#include <bits/stdc++.h>
#line 8 "library/gandalfr/other/io_supporter.hpp"
template <typename T>
std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
    for (int i = 0; i < (int)v.size(); i++)
        os << v[i] << (i + 1 != (int)v.size() ? " " : "");
    return os;
}
template <typename T>
std::ostream &operator<<(std::ostream &os, const std::set<T> &st) {
    for (const T &x : st) {
        std::cout << x << " ";
    }
    return os;
}
template <typename T>
std::ostream &operator<<(std::ostream &os, const std::multiset<T> &st) {
    for (const T &x : st) {
        std::cout << x << " ";
    }
    return os;
}
template <typename T>
std::ostream &operator<<(std::ostream &os, const std::deque<T> &dq) {
    for (const T &x : dq) {
        std::cout << x << " ";
    }
    return os;
}
template <typename T1, typename T2>
std::ostream &operator<<(std::ostream &os, const std::pair<T1, T2> &p) {
    os << p.first << ' ' << p.second;
    return os;
}
template <typename T>
std::ostream &operator<<(std::ostream &os, std::queue<T> &q) {
    int sz = q.size();
    while (--sz) {
        os << q.front() << ' ';
        q.push(q.front());
        q.pop();
    }
    os << q.front();
    q.push(q.front());
    q.pop();
    return os;
}
template <typename T>
std::istream &operator>>(std::istream &is, std::vector<T> &v) {
    for (T &in : v)
        is >> in;
    return is;
}
template <typename T1, typename T2>
std::istream &operator>>(std::istream &is, std::pair<T1, T2> &p) {
    is >> p.first >> p.second;
    return is;
}
#line 3 "library/gandalfr/math/matrix.hpp"
#line 8 "library/gandalfr/math/matrix.hpp"
template <class T> class matrix {
  private:
    int H, W;
    std::valarray<std::valarray<T>> table;
    enum rowtrans_operation_name { SCALE, SWAP, ADD };
    struct rowtrans_operation {
        int op, tar, res;
        T scl;
    };
    using operations_history = std::vector<rowtrans_operation>;
  public:
    matrix() = default;
    matrix(int _H, int _W, T val = 0)
        : H(_H), W(_W), table(std::valarray<T>(val, _W), _H) {}
    matrix(const std::vector<std::vector<T>> &vv)
        : H(vv.size()), W(vv[0].size()), table(std::valarray<T>(W), H) {
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                table[i][j] = vv[i][j];
    }
    matrix(const std::valarray<std::valarray<T>> &vv)
        : H(vv.size()), W(vv[0].size()), table(vv) {}
    /**
     * @brief 行列をリサイズする。
     * @param val 拡張部分の値
     */
    void resize(int _H, int _W, T val = 0) {
        H = _H, W = _W;
        table.resize(_H, std::valarray<T>(val, _H));
    }
    int size_H() const { return H; }
    int size_W() const { return W; }
    void transpose() {
        matrix<T> ret(W, H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                ret.table[j][i] = table[i][j];
        *this = std::move(ret);
    }
    void row_assign(int i, const std::valarray<T> &row) {
        assert(W == (int)row.size());
        table[i] = std::move(row);
    }
    void row_swap(int i, int j) {
        assert(0 <= i && i < H);
        assert(0 <= j && j < H);
        table[i].swap(table[j]);
    }
    /**
     * @attention O(n^3)
     * @attention 整数型では正しく計算できない。double や fraction を使うこと。
     * @attention 枢軸選びをしていないので double では誤差が出るかも。
     */
    operations_history sweep_method() {
        operations_history hist;
        T ret = 1;
        for (int h = 0, w = 0; h < H && w < W; w++) {
            if (table[h][w] == 0) {
                for (int piv = h + 1; piv < H; piv++) {
                    if (table[piv][w] != 0) {
                        hist.push_back({SWAP, h, piv, 0});
                        row_swap(h, piv);
                        break;
                    }
                }
                if (table[h][w] == 0) {
                    continue;
                }
            }
            T inv = 1 / table[h][w];
            hist.push_back({SCALE, -1, w, inv});
            table[h] *= inv;
            for (int j = h + 1; j < H; j++) {
                hist.push_back({ADD, h, j, -table[j][w]});
                table[j] -= table[h] * table[j][w];
            }
            h++;
        }
        return hist;
    }
    int rank() {
        auto U(*this);
        U.sweep_method();
        int r = 0;
        for (int i = 0; i < H; ++i) {
            for (int j = i; j < W; ++j) {
                if (U.table[i][j] != 0) {
                    r++;
                    break;
                }
            }
        }
        return r;
    }
    T determinant() const {
        assert(H == W);
        matrix<T> U(*this);
        T det = 1;
        auto hist = U.sweep_method();
        if (U.table[H - 1][H - 1] == 0)
            return 0;
        for (auto &[op, tar, res, scl] : hist) {
            switch (op) {
            case SCALE:
                det /= scl;
                break;
            case SWAP:
                det *= -1;
                break;
            }
        }
        return det;
    }
    std::vector<T> solve_system_of_equations(const std::vector<T> &y) {
        assert(H == W);
        std::vector<T> x(y);
        matrix<T> U(*this);
        auto hist = U.sweep_method();
        if (U.table[H - 1][H - 1] == 0)
            return {};
        for (auto &[op, tar, res, scl] : hist) {
            switch (op) {
            case SCALE:
                x[res] *= scl;
                break;
            case SWAP:
                std::swap(x[tar], x[res]);
                break;
            case ADD:
                x[res] += x[tar] * scl;
                break;
            }
        }
        for (int i = H - 1; i >= 0; --i) {
            for (int j = 0; j < i; ++j) {
                x[j] -= U.table[j][i] * x[i];
            }
        }
        return x;
    }
    matrix<T> inverse() {
        assert(H == W);
        matrix<T> INV(matrix<T>::E(H)), U(*this);
        auto hist = U.sweep_method();
        if (U.table[H - 1][H - 1] == 0)
            return matrix<T>(0, 0);
        for (auto &[op, tar, res, scl] : hist) {
            switch (op) {
            case SCALE:
                INV.table[res] *= scl;
                break;
            case SWAP:
                std::swap(INV.table[tar], INV.table[res]);
                break;
            case ADD:
                INV.table[res] += INV.table[tar] * scl;
                break;
            }
        }
        for (int i = H - 1; i >= 0; --i) {
            for (int j = 0; j < i; ++j) {
                INV.table[j] -= INV.table[i] * U.table[j][i];
            }
        }
        return INV;
    }
    void print() const {
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                std::cout << table[i][j] << (j == W - 1 ? "" : " ");
            }
            std::cout << std::endl;
        }
    }
    matrix<T> &operator+=(const matrix<T> &a) {
        this->table += a.table;
        return *this;
    }
    matrix<T> &operator-=(const matrix<T> &a) {
        this->table -= a.table;
        return *this;
    }
    matrix<T> &operator*=(const T &a) {
        this->table *= a;
        return *this;
    }
    matrix<T> &operator*=(const matrix<T> &a) {
        assert(W == a.H);
        matrix<T> a_t(a), ret(H, a.W);
        a_t.transpose();
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < a_t.H; j++) {
                ret.table[i][j] = (table[i] * a_t.table[j]).sum();
            }
        }
        *this = std::move(ret);
        return *this;
    }
    matrix<T> &operator/=(const T &a) {
        this->table /= a;
        return *this;
    }
    /**
     * @brief 行列の冪乗。
     * @param n 指数
     * @attention n が 0 なら単位行列。
     * @attention 演算子の優先度に注意。
     */
    matrix<T> operator^=(long long n) {
        assert(H == W);
        if (n == 0)
            return *this = E(H);
        n--;
        matrix<T> x(*this);
        while (n) {
            if (n & 1)
                *this *= x;
            x *= x;
            n >>= 1;
        }
        return *this;
    }
    matrix<T> operator+() { return *this; }
    matrix<T> operator-() { return matrix<T>(*this) *= -1; }
    matrix<T> operator+(const matrix<T> &a) { return matrix<T>(*this) += a; }
    matrix<T> operator-(const matrix<T> &a) { return matrix<T>(*this) -= a; }
    template <typename S> matrix<T> operator*(const S &a) {
        return matrix<T>(*this) *= a;
    }
    matrix<T> operator/(const T &a) { return matrix<T>(*this) /= a; }
    matrix<T> operator^(long long n) { return matrix<T>(*this) ^= n; }
    friend std::istream &operator>>(std::istream &is, matrix<T> &mt) {
        for (auto &arr : mt.table)
            for (auto &x : arr)
                is >> x;
        return is;
    }
    T &operator()(int h, int w) {
        assert(0 <= h && h < H && 0 <= w && w <= W);
        return table[h][w];
    }
    /**
     * @brief サイズ n の単位行列。
     */
    static matrix<T> E(int N) {
        matrix<T> ret(N, N);
        for (int i = 0; i < N; i++)
            ret.table[i][i] = 1;
        return ret;
    }
};
#line 4 "playspace/main.cpp"
using namespace std;
using ll = long long;
const int INF = 1001001001;
const int MAXINT = std::numeric_limits<int>::max();
const int MININT = std::numeric_limits<int>::min();
const ll INFLL = 1001001001001001001;
const ll MAXLL = std::numeric_limits<ll>::max();
const ll MINLL = std::numeric_limits<ll>::min();
const ll MOD  = 1000000007;
const ll _MOD = 998244353;
#define rep(i, j, n) for(ll i = (ll)(j); i < (ll)(n); i++)
#define rrep(i, j, n) for(ll i = (ll)(n-1); i >= (ll)(j); i--)
#define all(a) (a).begin(),(a).end()
#define debug(a) std::cerr << #a << ": " << a << std::endl
#define LF cout << endl
#ifdef ENABLE_MULTI_FOR
#define mrep(it, mr) for(multi_iter it(mr); !it.fin(); ++it)
#endif
template<typename T1, typename T2> inline bool chmax(T1 &a, T2 b) { return a < b && (a = b, true); }
template<typename T1, typename T2> inline bool chmin(T1 &a, T2 b) { return a > b && (a = b, true); }
void Yes(bool ok){ std::cout << (ok ? "Yes" : "No") << std::endl; }
int main(void) {
    matrix<ll> mt(2, 2);
    cin >> mt;
    (mt ^ 3).print();
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0