template <typename T = double, typename SFINAE = std::void_t<decltype(std::declval<T&>() * std::declval<T&>(), std::declval<T&>() + std::declval<T&>())>>
class matrix {
private:

    std::vector<std::vector<T>> data;

    size_t height;

    size_t width;

public:

    matrix(size_t n, size_t m) noexcept
            : data(m, std::vector<T>(n)),
              height(m),
              width(n)
              { }

    matrix(const matrix&) = default;
    matrix& operator=(const matrix&) = default;
    matrix(matrix&&) = default;
    matrix& operator=(matrix&&) = default;

    T& operator()(size_t x, size_t y) {
        if (x > width || y > height) {
            throw std::out_of_range("matrix index is too large");
        }
        return data[y][x];
    }

    const T operator()(size_t x, size_t y) const {
        if (x > width || y > height) {
            throw std::out_of_range("matrix index is too large");
        }
        return data[y][x];
    }

    constexpr size_t get_height() const {
        return height;
    }

    constexpr size_t get_width() const {
        return width;
    }

    matrix transpose() const {
        matrix result(width, height);
        for (size_t i = 0; i < height; ++i) {
            for (size_t j = 0; j < width; ++j) {
                result(i, j) = data[i][j];
            }
        }
        return result;
    }

};

template<typename T, typename U>
auto operator+(const matrix<T>& left, const matrix<U>& right) -> matrix<decltype(T() * U())> {
    if (left.get_width() != right.get_width() || left.get_height() != right.get_height()) {
        throw std::runtime_error("we need two matrix with the same dimension for addition");
    }
    matrix<decltype(T() * U())> result(left.get_height(), left.get_width());
    for (size_t i = 0; i < left.get_height(); ++i) {
        for (size_t j = 0; j < left.get_width(); ++j) {
            result(j, i) = left(j, i) + right(j, i);
        }
    }
    return result;
}

template<typename T, typename U>
auto operator*(const T multiplier, const matrix<U>& right) -> matrix<decltype(T() * U())> {
    matrix<decltype(T() * U())> result(right.get_height(), right.get_width());
    for (size_t i = 0; i < right.get_height(); ++i) {
        for (size_t j = 0; j < right.get_width(); ++j) {
            result(j, i) = multiplier * right(j, i);
        }
    }
    return result;
}

template<typename T, typename U>
auto operator*(const matrix<T>& left, const U multiplier) -> matrix<decltype(U() * T())> {
    matrix<decltype(T() * U())> result(left.get_height(), left.get_width());
    for (size_t i = 0; i < left.get_height(); ++i) {
        for (size_t j = 0; j < left.get_width(); ++j) {
            result(j, i) = left(j, i) * multiplier;
        }
    }
    return result;
}

template<typename T, typename U>
auto operator*(const matrix<T>& left, const matrix<U> right) -> matrix<decltype(U() * T())> {
    if (left.get_width() != right.get_height()) {
        throw std::runtime_error("cannot multiply two matrix if the dimensions aren't correct");
    }
    matrix<decltype(T() * U())> result(left.get_height(), right.get_width());
    for (size_t i = 0; i < left.get_height(); ++i) {
        for (size_t j = 0; j < right.get_width(); ++j) {
            for (size_t k = 0; k < left.get_width(); ++k) {
                result(j, i) += left(k, i) * right(j, k);
            }
        }
    }
    return result;
}

template<typename T, typename U>
auto hadamard_multiplicaton(const matrix<T>& left, const matrix<U>& right) -> matrix<decltype(std::declval<T&>() * std::declval<U&>())> {
    if (left.get_width() != right.get_width() || left.get_height() != right.get_height()) {
        throw std::runtime_error("we need two matrix with the same dimension for hadamard multiplication");
    }
    matrix<decltype(std::declval<T&>() * std::declval<U&>())> result(left.get_height(), left.get_width());
    for (size_t i = 0; i < left.get_height(); ++i) {
        for (size_t j = 0; j < left.get_width(); ++j) {
            result(j, i) = left(j, i) * right(j, i);
        }
    }
    return result;
}

template<typename T, typename U>
auto kronecker_product(const matrix<T>& left, const matrix<U>& right) -> matrix<decltype(std::declval<T&>() * std::declval<U&>())> {
    if (left.get_height() != 1 || right.get_width() != 1) {
        throw std::runtime_error("we need a row and a column vector for kronecker product");
    }
    matrix<decltype(std::declval<T&>() * std::declval<U&>())> result(right.get_height(), left.get_width());
    for (size_t i = 0; i < right.get_height(); ++i) {
        for (size_t j = 0; j < left.get_width(); ++j) {
            result(j, i) = left(j, 0) * right(0, i);
        }
    }
    return result;
}

template <typename T>
std::ostream& operator<<(std::ostream& os, const matrix<T>& m) {
    for (size_t i = 0; i < m.get_height(); ++i) {
        for (size_t j = 0; j < m.get_width(); ++j) {
            os << m(j, i) << "\t";
        }
        os << "\n";
    }
    return os;
}