gauss_jordan
(math/gauss_jordan.hpp)

• View this file on GitHub
• Last update: 2023-05-17 22:41:14+09:00
• Include: #include "math/gauss_jordan.hpp"

説明

vector<vector<T>> (Tは体) の行列に対して掃き出し法を適用したものを返す。 行列式などを得たい場合や、 double 型を扱う場合などには柔軟に実装する必要がある。

Depends on

• テンプレート (template/template.hpp)

Verified with

• test/others/aoj_1300.test.cpp

Code

#pragma once
#include "../template/template.hpp"

namespace lib {
using namespace std;

template <typename T> vector<vector<T>> gauss_jordan(vector<vector<T>> &a) {
    int n = a.size();
    int m = a[0].size();
    vector<vector<T>> b = a;
    int piv = 0;
    rep(j, 0, m) rep(i, piv, n) {
        if (b[i][j] != T(0)) {
            swap(b[i], b[piv]);
            T ip = T(1) / b[piv][j];
            rep(l, 0, n) {
                if (l != piv) {
                    T tmp = ip * b[l][j];
                    rep(k, j, m) b[l][k] -= tmp * b[piv][k];
                }
            }
            rep(k, j, m) b[piv][k] *= ip;
            piv++;
            break;
        }
    }
    return b;
}

}  // namespace lib

#line 2 "template/template.hpp"

#include <bits/stdc++.h>

#define rep(i, s, n) for (int i = (int)(s); i < (int)(n); i++)
#define rrep(i, s, n) for (int i = (int)(n)-1; i >= (int)(s); i--)
#define all(v) v.begin(), v.end()

using ll = long long;
using ld = long double;
using ull = unsigned long long;

template <typename T> bool chmin(T &a, const T &b) {
    if (a <= b) return false;
    a = b;
    return true;
}
template <typename T> bool chmax(T &a, const T &b) {
    if (a >= b) return false;
    a = b;
    return true;
}

namespace lib {

using namespace std;

}  // namespace lib

// using namespace lib;
#line 3 "math/gauss_jordan.hpp"

namespace lib {
using namespace std;

template <typename T> vector<vector<T>> gauss_jordan(vector<vector<T>> &a) {
    int n = a.size();
    int m = a[0].size();
    vector<vector<T>> b = a;
    int piv = 0;
    rep(j, 0, m) rep(i, piv, n) {
        if (b[i][j] != T(0)) {
            swap(b[i], b[piv]);
            T ip = T(1) / b[piv][j];
            rep(l, 0, n) {
                if (l != piv) {
                    T tmp = ip * b[l][j];
                    rep(k, j, m) b[l][k] -= tmp * b[piv][k];
                }
            }
            rep(k, j, m) b[piv][k] *= ip;
            piv++;
            break;
        }
    }
    return b;
}

}  // namespace lib

Back to top page