Learning Algorithms

アルゴリズムの勉強メモ

行列ライブラリ

行列ライブラリ

競プロ用に行列ライブラリを作った。

特別に綺麗だとか高速だとかそういうことはないけれど、もし利用したい人がいれば、勝手に利用していただいても構いません。多分動きます。

Typeは、値が小数になる場合と、逆行列が必要なときにdoubleにしておく。その他の場合はintでよい。
int GetRank(Matrix a)は任意の行列のランクを返す。
bool Inv(Matrix a, Matrix &inv)は第二引数に逆行列として受け取りたい、aと同じ大きさの空の行列を参照で渡して使う。もし逆行列が存在するならばtrueを返し、第二引数の中が逆行列に書き換えられる。存在しないならばfalseを返す。
AddSubは常に引数を2つ渡すようにすればよい。
いずれの関数も、定義できない演算をしようとすると、assertで怒られる。

また追記、修正するかもしれない。

typedef double Type;
typedef vector<vector<Type>> Matrix;
int GetRank(Matrix a) {
        int h = a.size(), w = a[0].size();
        int res = 0, now = 0;
        for (int i = 0; i < h; i ++) {
                Type ma = 0.0;
                int pivot;
                for (int j = i; j < h; j ++) {
                        if (a[j][now] > ma) {
                                ma = a[j][now];
                                pivot = j;
                        }
                }
                if (ma == 0.0) {
                        now ++;
                        if (now == w) break;
                        i --;
                        continue;
                }
                if (pivot != i) {
                        for (int j = 0; j < w; j ++) { 
                                swap(a[i][j], a[pivot][j]);
                        }

                }
                Type tmp = 1.0 / a[i][now];
                for (int j = 0; j < w; j ++) a[i][j] *= tmp;
                for (int j = 0; j < h; j ++) {
                        if (i != j) {
                                Type tmp2 = a[j][now];
                                for (int k = 0; k < w; k ++) {
                                        a[j][k] -= a[i][k] * tmp2;
                                }
                        }
                }
                res ++;
        }
        return res;
}
bool Inv(Matrix a, Matrix &inv) {
        assert(a.size() == a[0].size() && inv.size() == inv[0].size());
        int n = a.size();
        for (int i = 0; i < n; i ++) {
                for (int j = 0; j < n; j ++) {
                        inv[i][j] = (i == j ? 1.0 : 0.0);
                }
        }
        for (int i = 0; i < n; i ++) {
                Type ma = 0.0;
                int pivot;
                for (int j = i; j < n; j ++) { 
                        if (a[j][i] > ma) {
                                ma = a[j][i];
                                pivot = j;
                        }
                }
                if (ma == 0.0) return false;
                if (pivot != i) {
                        for (int j = 0; j < n; j ++) { 
                                swap(a[i][j], a[pivot][j]);
                                swap(inv[i][j], inv[pivot][j]);
                        }

                }
                Type tmp = 1.0 / a[i][i];
                for (int j = 0; j < n; j ++) {
                        a[i][j] *= tmp;
                        inv[i][j] *= tmp;
                }
                for (int j = 0; j < n; j ++) {
                        if (i != j) {
                                Type tmp2 = a[j][i];
                                for (int k = 0; k < n; k ++) {
                                        a[j][k] -= a[i][k] * tmp2;
                                        inv[j][k] -= inv[i][k] * tmp2;
                                }
                        }
                }
        }
        return true;
}
Matrix Add(const Matrix &a, const Matrix &b, bool minus = false) {
        assert(a.size() == b.size() && a[0].size() == b[0].size());
        int h = a.size(), w = a[0].size();
        Matrix c(h, vector<Type> (w));
        for (int i = 0; i < h; i ++) {
                for (int j = 0; j < w; j ++) {
                        c[i][j] = a[i][j] + (minus ? -1 : 1) * b[i][j];
                }
        }
        return c;
}
Matrix Sub(const Matrix &a, const Matrix &b) {
        return Add(a, b, true);
}
Matrix Mul(const Matrix &a, const Matrix &b) {
        assert(a[0].size() == b.size());
        Matrix c(a.size(), vector<Type> (b[0].size()));
        for (int i = 0; i < a.size(); i ++) {
                for (int k = 0; k < b.size(); k ++) {
                        for (int j = 0; j < b[0].size(); j ++) {
                                c[i][j] = (c[i][j] + a[i][k] * b[k][j]);
                        }
                }
        }
        return c;
}
Matrix Pow(Matrix a, long long n) {
        assert(a.size() == a[0].size());
        Matrix b(a.size(), vector<Type> (a.size()));
        for (int i = 0; i < a.size(); i ++) {
                b[i][i] = 1;
        }
        while (n > 0) {
                if (n & 1) b = Mul(b, a);
                a = Mul(a, a);
                n >>= 1;
        }
        return b;
}
void PrintMatrix(const Matrix &a) {
        int h = a.size(), w = a[0].size();
        for (int i = 0; i < h; i ++) {
                for (int j = 0; j < w; j ++) {
                        cout << a[i][j] << ' ';
                }
                cout << endl;
        }
}