diff --git a/c++/nda/mapped_functions.hpp b/c++/nda/mapped_functions.hpp
index ec43a1b86..cb9b22812 100644
--- a/c++/nda/mapped_functions.hpp
+++ b/c++/nda/mapped_functions.hpp
@@ -67,7 +67,7 @@ namespace nda {
   /**
    * @brief Function pow for nda::ArrayOrScalar types (lazy and coefficient-wise for nda::Array types).
    *
-   * @tparam A nda::ArrayOrScalar type..
+   * @tparam A nda::ArrayOrScalar type.
    * @param a nda::ArrayOrScalar object.
    * @param p Exponent value.
    * @return A lazy nda::expr_call object (nda::Array) or the result of `std::pow` applied to the object (nda::Scalar).
@@ -84,7 +84,7 @@ namespace nda {
    * @brief Function conj for nda::ArrayOrScalar types (lazy and coefficient-wise for nda::Array types with a complex
    * value type).
    *
-   * @tparam A nda::ArrayOrScalar type..
+   * @tparam A nda::ArrayOrScalar type.
    * @param a nda::ArrayOrScalar object.
    * @return A lazy nda::expr_call object (nda::Array and complex valued), the forwarded input object (nda::Array and
    * not complex valued) or the complex conjugate of the scalar input.
@@ -97,6 +97,21 @@ namespace nda {
       return std::forward<A>(a);
   }
 
+  /**
+   * @brief Reciprocal function for nda::ArrayOrScalar types (lazy and coefficient-wise for nda::Array types).
+   *
+   * @tparam A nda::ArrayOrScalar type.
+   * @param a nda::ArrayOrScalar object.
+   * @return A lazy nda::expr_call object (nda::Array) or the result of \f$ 1 / x \f$ applied to the object 
+   * (nda::Scalar).
+   */
+  template <ArrayOrScalar A>
+  auto reciprocal(A &&a) {
+    return nda::map([](auto const &x) {
+      return 1 / x;
+    })(std::forward<A>(a));
+  }
+
   /** @} */
 
 } // namespace nda
diff --git a/test/c++/nda_math_functions.cpp b/test/c++/nda_math_functions.cpp
index 1a6b12d0d..51a1c40a7 100644
--- a/test/c++/nda_math_functions.cpp
+++ b/test/c++/nda_math_functions.cpp
@@ -228,6 +228,17 @@ TEST_F(NDAMathFunction, Real) {
   EXPECT_EQ(nda::real(cplx_c), std::real(cplx_c));
 }
 
+TEST_F(NDAMathFunction, Reciprocal) {
+  auto B_d = nda::reciprocal(A_d);
+  auto B_c = nda::reciprocal(A_c);
+  nda::for_each(shape, [&](auto... idxs) {
+    EXPECT_DOUBLE_EQ(B_d(idxs...), 1 / A_d(idxs...));
+    EXPECT_COMPLEX_NEAR(B_c(idxs...), 1 / A_c(idxs...), 1e-14);
+  });
+  EXPECT_DOUBLE_EQ(nda::reciprocal(dbl_c), 1 / dbl_c);
+  EXPECT_COMPLEX_NEAR(nda::reciprocal(cplx_c), 1 / cplx_c);
+}
+
 TEST_F(NDAMathFunction, Sin) {
   auto B_d = nda::sin(A_d);
   auto B_c = nda::sin(A_c);
@@ -295,10 +306,12 @@ TEST_F(NDAMathFunction, Combinations) {
   auto B_d = nda::pow(nda::pow(nda::abs(A_d), 1.5), 2.0 / 3.0);
   auto C_d = nda::sqrt(nda::pow(A_d, 2));
   auto D_d = nda::log(nda::exp(A_d));
+  auto B_c = nda::reciprocal(nda::reciprocal(A_c));
   nda::for_each(shape, [&](auto... idxs) {
     EXPECT_NEAR(B_d(idxs...), std::abs(A_d(idxs...)), 1e-10);
     EXPECT_NEAR(C_d(idxs...), std::abs(A_d(idxs...)), 1e-10);
     EXPECT_NEAR(D_d(idxs...), A_d(idxs...), 1e-10);
+    EXPECT_COMPLEX_NEAR(B_c(idxs...), A_c(idxs...), 1e-10);
   });
 }