#include <iostream>
#include <vector>

class Matrix {
public:
    Matrix(size_t size) : _size(size), _data(size * size) {}
    Matrix(size_t size, std::vector<double> data) : _size(size), _data(size * size) {
        _data = data;
    }
    double& operator()(size_t i, size_t j) {
        return _data[i * _size + j];
    }
    double operator()(size_t i, size_t j) const {
        return _data[i * _size + j];
    }
    Matrix I(size_t size) {
        auto matrix = Matrix(size);
        for (size_t i = 0; i < size; ++i)
            matrix(i, i) = 1;

        return matrix;
    }

    Matrix operator+(const Matrix& b) {
        auto sum = Matrix(_size, _data);
        for (size_t i = 0; i < _size; ++i)
            for (size_t j = 0; j < _size; ++j)
                sum(i, j) += b(i, j);

        return sum;
    }

    Matrix& operator+=(const Matrix& b) {
        for (size_t i = 0; i < _size; ++i)
            for (size_t j = 0; j < _size; ++j)
                operator()(i, j) += b(i, j);

        return *this;
    }

    Matrix operator*(const Matrix& b) {
        auto product = Matrix(_size);

        for (size_t i = 0; i < _size; ++i)
            for (size_t j = 0; j < _size; ++j)
                for (size_t k = 0; k < _size; ++k) {
                    product(i, j) += operator()(i, k) * b(k, j);
                }

        return product;
    }

    Matrix operator*(const double& d) {
        auto scalarMult = Matrix(_size, _data);
        for (size_t i = 0; i < _size; ++i)
            for (size_t j = 0; j < _size; ++j)
                scalarMult(i, j) *= d;

        return scalarMult;
    }

    Matrix pow(size_t n) {
        auto result = Matrix(_size);

        for (size_t i = 0; i < _size; i++)
        {
            result(i, i) = 1;
        }

        for (size_t i = 0; i < n; i++) {
            result = result * *this;
        }

        return result;
    }


private:
    size_t _size;
    std::vector<double> _data;
};

int main() {
    auto A = Matrix(2, {1, 1, 1, 0});
    auto P1 = Matrix(2);
    auto P2 = Matrix(2, {1, 1, 1, 0});

    std::vector<double> polinomialCoefs { 1, 2, 3, 4, 5 };

    auto test = A.pow(0);

    for (int i = 0; i < 5; i++) {
        P1 += A.pow(i) * polinomialCoefs[i];
    }

    P2 = Matrix.I(2);
    for (int i = 4; i >= 1; i--) {
        P2 += A * P2 + polinomialCoefs[i - 1];
    }

    return 0;
}