#include <iostream>
#include <vector>
#include <stdexcept>
#include <initializer_list>
#include <type_traits>

// Template class for Matrix
template <typename T>
class Matrix {
public:
    Matrix(size_t rows, size_t cols)
        : rows_(rows), cols_(cols), data_(rows, std::vector<T>(cols)) {}

    Matrix(std::initializer_list<std::initializer_list<T>> init) {
        rows_ = init.size();
        cols_ = init.begin()->size();
        data_.resize(rows_);
        size_t i = 0;
        for (const auto& row : init) {
            data_[i] = row;
            ++i;
        }
    }

    size_t rows() const { return rows_; }
    size_t cols() const { return cols_; }

    T& operator()(size_t row, size_t col) {
        if (row >= rows_ || col >= cols_) {
            throw std::out_of_range("Matrix indices out of range");
        }
        return data_[row][col];
    }

    const T& operator()(size_t row, size_t col) const {
        return const_cast<Matrix&>(*this)(row, col);
    }

    Matrix<T> operator*(const Matrix<T>& other) const {
        if (cols_ != other.rows()) {
            throw std::invalid_argument("Matrix dimensions do not match for multiplication");
        }
        Matrix<T> result(rows_, other.cols());
        for (size_t i = 0; i < rows_; ++i) {
            for (size_t j = 0; j < other.cols(); ++j) {
                for (size_t k = 0; k < cols_; ++k) {
                    result(i, j) += (*this)(i, k) * other(k, j);
                }
            }
        }
        return result;
    }

    void print() const {
        for (const auto& row : data_) {
            for (const auto& elem : row) {
                std::cout << elem << " ";
            }
            std::cout << std::endl;
        }
    }

private:
    size_t rows_, cols_;
    std::vector<std::vector<T>> data_;
};

// Variadic template function for matrix multiplication
template<typename T>
Matrix<T> matrixProduct(const Matrix<T>& m) {
    return m;
}

template<typename T, typename... Args>
Matrix<T> matrixProduct(const Matrix<T>& first, const Args&... args) {
    return first * matrixProduct(args...);
}

// Constexpr for physical constants
constexpr double SPEED_OF_LIGHT = 299792458.0; // m/s
constexpr double GRAVITATIONAL_CONSTANT = 6.67430e-11; // N(m/kg)^2

// Function to demonstrate decltype usage
template<typename T, typename U>
auto addMatrices(const Matrix<T>& m1, const Matrix<U>& m2) -> Matrix<decltype(std::declval<T>() + std::declval<U>())> {
    if (m1.rows() != m2.rows() || m1.cols() != m2.cols()) {
        throw std::invalid_argument("Matrix dimensions do not match for addition");
    }
    Matrix<decltype(std::declval<T>() + std::declval<U>())> result(m1.rows(), m1.cols());
    for (size_t i = 0; i < m1.rows(); ++i) {
        for (size_t j = 0; j < m1.cols(); ++j) {
            result(i, j) = m1(i, j) + m2(i, j);
        }
    }
    return result;
}

int main() {
    // Demonstrate Matrix template class
    Matrix<int> m1{{1, 2}, {3, 4}};
    Matrix<int> m2{{5, 6}, {7, 8}};

    std::cout << "Matrix m1:" << std::endl;
    m1.print();
    std::cout << "Matrix m2:" << std::endl;
    m2.print();

    // Demonstrate variadic template matrix multiplication
    auto result = matrixProduct(m1, m2, m1);
    std::cout << "Result of m1 * m2 * m1:" << std::endl;
    result.print();

    // Demonstrate constexpr usage
    std::cout << "Speed of light: " << SPEED_OF_LIGHT << " m/s" << std::endl;
    std::cout << "Gravitational constant: " << GRAVITATIONAL_CONSTANT << " N(m/kg)^2" << std::endl;

    // Demonstrate decltype usage
    Matrix<double> m3{{1.5, 2.5}, {3.5, 4.5}};
    auto sum = addMatrices(m1, m3);
    std::cout << "Sum of integer and double matrices:" << std::endl;
    sum.print();

    return 0;
}