From e116f6847e3cd9d16ea3be14c2b7efbd0a1c2b0b Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 4 Apr 2018 13:47:23 +0200
Subject: bug #1521: avoid signalling NaN in hypot and make it std::complex<>
 friendly.

---
 Eigen/src/Core/MathFunctions.h     | 26 +-------------------------
 Eigen/src/Core/MathFunctionsImpl.h | 19 +++++++++++++++++++
 2 files changed, 20 insertions(+), 25 deletions(-)

(limited to 'Eigen')

diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index e981129b2..84f9d0cd5 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -347,31 +347,7 @@ struct norm1_retval
 * Implementation of hypot                                                *
 ****************************************************************************/
 
-template<typename Scalar>
-struct hypot_impl
-{
-  typedef typename NumTraits<Scalar>::Real RealScalar;
-  static inline RealScalar run(const Scalar& x, const Scalar& y)
-  {
-    EIGEN_USING_STD_MATH(abs);
-    EIGEN_USING_STD_MATH(sqrt);
-    RealScalar _x = abs(x);
-    RealScalar _y = abs(y);
-    Scalar p, qp;
-    if(_x>_y)
-    {
-      p = _x;
-      qp = _y / p;
-    }
-    else
-    {
-      p = _y;
-      qp = _x / p;
-    }
-    if(p==RealScalar(0)) return RealScalar(0);
-    return p * sqrt(RealScalar(1) + qp*qp);
-  }
-};
+template<typename Scalar> struct hypot_impl;
 
 template<typename Scalar>
 struct hypot_retval
diff --git a/Eigen/src/Core/MathFunctionsImpl.h b/Eigen/src/Core/MathFunctionsImpl.h
index ae1386b4c..034cfad7b 100644
--- a/Eigen/src/Core/MathFunctionsImpl.h
+++ b/Eigen/src/Core/MathFunctionsImpl.h
@@ -66,6 +66,25 @@ T generic_fast_tanh_float(const T& a_x)
   return pdiv(p, q);
 }
 
+
+template<typename Scalar>
+struct hypot_impl
+{
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  static inline RealScalar run(const Scalar& x, const Scalar& y)
+  {
+    EIGEN_USING_STD_MATH(abs);
+    EIGEN_USING_STD_MATH(sqrt);
+    RealScalar _x = abs(x);
+    RealScalar _y = abs(y);
+    RealScalar p, qp;
+    p = numext::maxi(_x,_y);
+    if(p==RealScalar(0)) return RealScalar(0);
+    qp = numext::mini(_y,_x) / p;    
+    return p * sqrt(RealScalar(1) + qp*qp);
+  }
+};
+
 } // end namespace internal
 
 } // end namespace Eigen
-- 
cgit v1.2.3