Number/fixed_matrix.hpp

• View this file on GitHub
• Last update: 2020-11-05 12:15:41+09:00
• Include: #include "Number/fixed_matrix.hpp"

Verified with

• Verify/fixed_matrix_pow.test.cpp

Code

#pragma once

#include <array>

template <class T, int D>
struct Matrix : public std::array<std::array<T, D>, D> {
    // constructor
    using std::array<std::array<T, D>, D>::array;

    Matrix(T val = 0) {
        for (auto& v : *this) v.fill(val);
    }

    static Matrix id() {
        Matrix m;
        for (int i = 0; i < D; ++i) m[i][i] = 1;
        return m;
    }

    // arithmetic
    Matrix operator*(const Matrix& m) const { return Matrix(*this) *= m; }
    Matrix& operator*=(const Matrix& m) {
        Matrix nmat;
        for (auto& v : nmat) v.fill(0);

        for (int i = 0; i < D; ++i) {
            for (int j = 0; j < D; ++j) {
                for (int k = 0; k < D; ++k) {
                    nmat[i][j] += (*this)[i][k] * m[k][j];
                }
            }
        }
        return *this = nmat;
    }

    template <class U>
    Matrix pow(U k) {
        Matrix ret = id();
        Matrix a = *this;

        while (k > 0) {
            if (k & 1) ret *= a;
            a *= a;
            k >>= 1;
        }
        return ret;
    }
};

#line 2 "Number/fixed_matrix.hpp"

#include <array>

template <class T, int D>
struct Matrix : public std::array<std::array<T, D>, D> {
    // constructor
    using std::array<std::array<T, D>, D>::array;

    Matrix(T val = 0) {
        for (auto& v : *this) v.fill(val);
    }

    static Matrix id() {
        Matrix m;
        for (int i = 0; i < D; ++i) m[i][i] = 1;
        return m;
    }

    // arithmetic
    Matrix operator*(const Matrix& m) const { return Matrix(*this) *= m; }
    Matrix& operator*=(const Matrix& m) {
        Matrix nmat;
        for (auto& v : nmat) v.fill(0);

        for (int i = 0; i < D; ++i) {
            for (int j = 0; j < D; ++j) {
                for (int k = 0; k < D; ++k) {
                    nmat[i][j] += (*this)[i][k] * m[k][j];
                }
            }
        }
        return *this = nmat;
    }

    template <class U>
    Matrix pow(U k) {
        Matrix ret = id();
        Matrix a = *this;

        while (k > 0) {
            if (k & 1) ret *= a;
            a *= a;
            k >>= 1;
        }
        return ret;
    }
};

Back to top page