From 8f1b6198c26d8c72bdb7762ac56611eb06354f30 Mon Sep 17 00:00:00 2001
From: Chip Kerchner <chip.kerchner@ibm.com>
Date: Thu, 16 Feb 2023 23:35:42 +0000
Subject: [PATCH] Fix epsilon and dummy_precision values in long double for
 double doubles.  Prevented some algorithms from converging on PPC.

(cherry picked from commit 54459214a1b9c67df04bc529474fca1ec9f4c84f)
---
 Eigen/src/Core/NumTraits.h                         | 14 ++++++++++++--
 .../Eigen/src/MatrixFunctions/MatrixPower.h        |  2 +-
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h
index 72eac5a93..8e81237ac 100644
--- a/Eigen/src/Core/NumTraits.h
+++ b/Eigen/src/Core/NumTraits.h
@@ -248,8 +248,18 @@ template<> struct NumTraits<double> : GenericNumTraits<double>
 template<> struct NumTraits<long double>
   : GenericNumTraits<long double>
 {
-  EIGEN_CONSTEXPR
-  static inline long double dummy_precision() { return 1e-15l; }
+  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
+  static inline long double dummy_precision() { return static_cast<long double>(1e-15l); }
+
+#if defined(EIGEN_ARCH_PPC) && (__LDBL_MANT_DIG__ == 106)
+  // PowerPC double double causes issues with some values
+  EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
+  static inline long double epsilon()
+  {
+    // 2^(-(__LDBL_MANT_DIG__)+1)
+    return static_cast<long double>(2.4651903288156618919116517665087e-32l);
+  }
+#endif
 };
 
 template<typename _Real> struct NumTraits<std::complex<_Real> >
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index d7672d7c9..ce92f5bfd 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -281,7 +281,7 @@ inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT)
 #endif
   int degree = 3;
   for (; degree <= maxPadeDegree; ++degree)
-    if (normIminusT <= maxNormForPade[degree - 3])
+    if (normIminusT <= static_cast<long double>(maxNormForPade[degree - 3]))
       break;
   return degree;
 }