From e25e3a041bc41266df31f48a441d95b38fc5a6fc Mon Sep 17 00:00:00 2001
From: Benoit Steiner <benoit.steiner.goog@gmail.com>
Date: Thu, 3 Dec 2015 18:16:35 -0800
Subject: [PATCH] Added rsqrt() method to the Array class: this method computes
 the coefficient-wise inverse square root much more efficiently than calling
 sqrt().inverse().

---
 Eigen/src/plugins/ArrayCwiseUnaryOps.h | 17 +++++++++++++++++
 test/array.cpp                         |  5 +++++
 2 files changed, 22 insertions(+)

diff --git a/Eigen/src/plugins/ArrayCwiseUnaryOps.h b/Eigen/src/plugins/ArrayCwiseUnaryOps.h
index a9310f12d..45e826b0c 100644
--- a/Eigen/src/plugins/ArrayCwiseUnaryOps.h
+++ b/Eigen/src/plugins/ArrayCwiseUnaryOps.h
@@ -4,6 +4,7 @@ typedef CwiseUnaryOp<internal::scalar_abs_op<Scalar>, const Derived> AbsReturnTy
 typedef CwiseUnaryOp<internal::scalar_arg_op<Scalar>, const Derived> ArgReturnType;
 typedef CwiseUnaryOp<internal::scalar_abs2_op<Scalar>, const Derived> Abs2ReturnType;
 typedef CwiseUnaryOp<internal::scalar_sqrt_op<Scalar>, const Derived> SqrtReturnType;
+typedef CwiseUnaryOp<internal::scalar_rsqrt_op<Scalar>, const Derived> RsqrtReturnType;
 typedef CwiseUnaryOp<internal::scalar_sign_op<Scalar>, const Derived> SignReturnType;
 typedef CwiseUnaryOp<internal::scalar_inverse_op<Scalar>, const Derived> InverseReturnType;
 typedef CwiseUnaryOp<internal::scalar_boolean_not_op<Scalar>, const Derived> BooleanNotReturnType;
@@ -139,6 +140,22 @@ sqrt() const
   return SqrtReturnType(derived());
 }
 
+/** \returns an expression of the coefficient-wise inverse square root of *this.
+  *
+  * This function computes the coefficient-wise inverse square root.
+  *
+  * Example: \include Cwise_sqrt.cpp
+  * Output: \verbinclude Cwise_sqrt.out
+  *
+  * \sa pow(), square()
+  */
+EIGEN_DEVICE_FUNC
+inline const RsqrtReturnType
+rsqrt() const
+{
+  return RsqrtReturnType(derived());
+}
+
 /** \returns an expression of the coefficient-wise signum of *this.
   *
   * This function computes the coefficient-wise signum.
diff --git a/test/array.cpp b/test/array.cpp
index 367bda2c4..5395721f5 100644
--- a/test/array.cpp
+++ b/test/array.cpp
@@ -236,6 +236,7 @@ template<typename ArrayType> void array_real(const ArrayType& m)
   // avoid NaNs with abs() so verification doesn't fail
   m3 = m1.abs();
   VERIFY_IS_APPROX(m3.sqrt(), sqrt(abs(m1)));
+  VERIFY_IS_APPROX(m3.rsqrt(), Scalar(1)/sqrt(abs(m1)));
   VERIFY_IS_APPROX(m3.log(), log(m3));
   VERIFY_IS_APPROX(m3.log10(), log10(m3));
 
@@ -292,6 +293,10 @@ template<typename ArrayType> void array_real(const ArrayType& m)
 
   VERIFY_IS_APPROX(m3.pow(RealScalar(0.5)), m3.sqrt());
   VERIFY_IS_APPROX(pow(m3,RealScalar(0.5)), m3.sqrt());
+
+  VERIFY_IS_APPROX(m3.pow(RealScalar(-0.5)), m3.rsqrt());
+  VERIFY_IS_APPROX(pow(m3,RealScalar(-0.5)), m3.rsqrt());
+
   VERIFY_IS_APPROX(log10(m3), log(m3)/log(10));
 
   // scalar by array division