cpp23.cpp

// C++23 Linear Algebra Example: R = A * B + c
#include <array>
#include <iostream>
#include <iomanip>

template<size_t Rows, size_t Cols>
class Matrix {
public:
    std::array<std::array<double, Cols>, Rows> data{};

    constexpr double& operator()(size_t i, size_t j) noexcept {
        return data[i][j];
    }
    constexpr const double& operator()(size_t i, size_t j) const noexcept {
        return data[i][j];
    }
};

// Matrix multiplication
template<size_t M, size_t K, size_t N> constexpr auto operator*(
    const Matrix<M, K>& A, const Matrix<K, N>& B
) {
    Matrix<M, N> R{};
    for (size_t i = 0; i < M; ++i)
        for (size_t j = 0; j < N; ++j)
            for (size_t k = 0; k < K; ++k)
                R(i, j) += A(i, k) * B(k, j);

    return R;
}

// Scalar addition
template<size_t Rows, size_t Cols> constexpr auto operator+(
    const Matrix<Rows, Cols>& matrix, double val
) {
    auto result = matrix;
    for (size_t i = 0; i < Rows; ++i)
        for (size_t j = 0; j < Cols; ++j)
            result(i, j) += val;

    return result;
}

// Print matrix
template<size_t Rows, size_t Cols>
std::ostream& operator<<(
    std::ostream& os, const Matrix<Rows, Cols>& matrix
) {
    for (size_t i = 0; i < Rows; ++i) {
        for (size_t j = 0; j < Cols; ++j) {
            os << std::fixed << std::setprecision(2) << matrix(i, j) << "\t";
        }

        os << '\n';
    }
    return os;
}

int main() {
    constexpr std::size_t M_DIM = 2;
    constexpr std::size_t K_DIM = 2;
    constexpr std::size_t N_DIM = 4;

    constexpr Matrix<M_DIM, K_DIM> A{{
        std::array{1.0, 2.0},
        std::array{3.0, 4.0}
    }};

    constexpr Matrix<K_DIM, N_DIM> B{{
        std::array{5.0, 6.0, 7.0, 8.0},
        std::array{9.0, 10.0, 11.0, 12.0}
    }};

    constexpr double c = 5.0;

    constexpr auto R = A * B + c;

    std::cout << "R = A * B + c:\n" << R;
}