diff --git a/include/rotgen/impl/matrix_impl64.hpp b/include/rotgen/impl/matrix_impl64.hpp
index 0764624..1651f71 100644
--- a/include/rotgen/impl/matrix_impl64.hpp
+++ b/include/rotgen/impl/matrix_impl64.hpp
@@ -43,6 +43,15 @@ namespace rotgen
     void transposeInPlace();
     void adjointInPlace();
 
+    double sum() const;
+    double prod() const;
+    double mean() const;
+    double maxCoeff() const;
+    double maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
+    double minCoeff() const;
+    double minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const;
+    double trace() const;
+
     double&       operator()(std::size_t i, std::size_t j);
     double const& operator()(std::size_t i, std::size_t j) const;
 
diff --git a/src/matrix_impl64.cpp b/src/matrix_impl64.cpp
index 40d616b..34cb985 100644
--- a/src/matrix_impl64.cpp
+++ b/src/matrix_impl64.cpp
@@ -101,11 +101,24 @@ namespace rotgen
   {
     storage_->data.transposeInPlace();
   }
+
   void matrix_impl64::adjointInPlace()
   {
     storage_->data.adjointInPlace();
   }
 
+  double matrix_impl64::sum() const { return storage_->data.sum(); }
+  double matrix_impl64::prod() const { return storage_->data.prod(); }
+  double matrix_impl64::mean() const { return storage_->data.mean(); }
+
+  double matrix_impl64::minCoeff() const { return storage_->data.minCoeff(); }
+  double matrix_impl64::minCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.minCoeff(row, col); }
+
+  double matrix_impl64::maxCoeff() const { return storage_->data.maxCoeff(); }
+  double matrix_impl64::maxCoeff(std::ptrdiff_t* row, std::ptrdiff_t* col) const { return storage_->data.maxCoeff(row, col); }
+
+  double matrix_impl64::trace() const { return storage_->data.trace(); }
+
   //==================================================================================================
   // Operators
   //==================================================================================================
diff --git a/test/basic/arithmetic_functions.cpp b/test/basic/arithmetic_functions.cpp
index 8d8ad26..c61f505 100644
--- a/test/basic/arithmetic_functions.cpp
+++ b/test/basic/arithmetic_functions.cpp
@@ -8,6 +8,7 @@
 #define TTS_MAIN
 #include <rotgen/matrix.hpp>
 #include "tts.hpp"
+#include <Eigen/Dense>
 
 template <typename MatrixType>
 struct MatrixDescriptor
@@ -44,6 +45,46 @@ void test_matrix_unary_ops(std::size_t rows, std::size_t cols,
   TTS_EQUAL(original, transposed_matrix);
 }
 
+template <typename MatrixType>
+void test_matrix_scalar_reductions(std::size_t rows, std::size_t cols,
+                                   const std::function<void(MatrixType&, std::size_t, std::size_t)>& init_fn)
+{
+  using EigenMatrix = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
+
+  MatrixType original(rows, cols);
+  EigenMatrix ref(rows, cols);
+
+  for (std::size_t r = 0; r < rows; ++r)
+    for (std::size_t c = 0; c < cols; ++c)
+    {
+      init_fn(original, r, c);
+      ref(r, c) = original(r, c);
+    }
+
+  TTS_EQUAL(original.sum(), ref.sum());
+  TTS_EQUAL(original.prod(), ref.prod());
+  TTS_EQUAL(original.mean(), ref.mean());
+  TTS_EQUAL(original.maxCoeff(), ref.maxCoeff());
+  TTS_EQUAL(original.minCoeff(), ref.minCoeff());
+  TTS_EQUAL(original.trace(), ref.trace());
+
+  std::ptrdiff_t row, col, ref_row, ref_col;
+
+  double cmin = original.minCoeff(&row, &col);
+  double rmin = ref.minCoeff(&ref_row, &ref_col);
+
+  TTS_EQUAL(cmin, rmin);
+  TTS_EQUAL(row, ref_row);
+  TTS_EQUAL(col, ref_col);
+
+  double cmax = original.maxCoeff(&row, &col);
+  double rmax = ref.maxCoeff(&ref_row, &ref_col);
+
+  TTS_EQUAL(cmax, rmax);
+  TTS_EQUAL(row, ref_row);
+  TTS_EQUAL(col, ref_col);
+}
+
 TTS_CASE("Matrix unary operations: transpose, adjoint, conjugate")
 {
   std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
@@ -61,3 +102,20 @@ TTS_CASE("Matrix unary operations: transpose, adjoint, conjugate")
     test_matrix_unary_ops<rotgen::matrix<double>>(desc.rows, desc.cols, desc.init_fn);
 };
 
+TTS_CASE("Basic arithmetic reduction operations")
+{
+  std::vector<MatrixDescriptor<rotgen::matrix<double>>> test_matrices = {
+    {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = r + c; }},
+    {3, 3, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 0.0; }},
+    {2, 4, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -r -c*c - 1234; }},
+    {4, 4, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 7.0; }},
+    {1, 1, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = 42.0; }},
+    {4, 2, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = std::sin(r + c); }},
+    {1, 5, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = -1.5 * r + 2.56 * c; }},
+    {5, 7, [](auto& m, std::size_t r, std::size_t c) { m(r, c) = (r == c ? 1.0 : 0.0); }},
+  };
+
+  for (const auto& [rows, cols, init_fn] : test_matrices)
+    test_matrix_scalar_reductions<rotgen::matrix<double>>(rows, cols, init_fn);
+};
+