From 95c0769e29c28ca414bf4644f2ce73b47a32ae5e Mon Sep 17 00:00:00 2001
From: Christopher Kormanyos <ckormanyos@yahoo.com>
Date: Mon, 1 Feb 2021 09:15:21 +0100
Subject: [PATCH] [CI SKIP] Handle more special cases and test zero

---
 .../math/cstdfloat/cstdfloat_complex_std.hpp  | 76 +++++++++++++++++--
 test/float128/test_std_lib.cpp                |  4 +
 2 files changed, 72 insertions(+), 8 deletions(-)

diff --git a/include/boost/math/cstdfloat/cstdfloat_complex_std.hpp b/include/boost/math/cstdfloat/cstdfloat_complex_std.hpp
index 1dd9ebc79..697999f04 100644
--- a/include/boost/math/cstdfloat/cstdfloat_complex_std.hpp
+++ b/include/boost/math/cstdfloat/cstdfloat_complex_std.hpp
@@ -509,13 +509,61 @@
     inline complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> pow(const BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE& x,
                                                                     const complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>& a)
     {
+      const bool x_isneg = (x < 0);
+      const bool x_isnan = (x.real() != x.real());
+      const bool x_isinf = ((!x_isneg) ? bool(+x.real() > (std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::max)())
+                                        : bool(-x.real() > (std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::max)()));
+
+      const bool a_re_isneg = (a.real() < 0);
+      const bool a_re_isnan = (a.real() != a.real());
+      const bool a_re_isinf = ((!a_re_isneg) ? bool(+a.real() > (std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::max)())
+                                             : bool(-a.real() > (std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::max)()));
+
+      const bool a_im_isneg = (a.imag() < 0);
+      const bool a_im_isnan = (a.imag() != a.imag());
+      const bool a_im_isinf = ((!a_im_isneg) ? bool(+a.imag() > (std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::max)())
+                                             : bool(-a.imag() > (std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::max)()));
+
+      const bool args_is_nan = (x_isnan || a_re_isnan || a_im_isnan);
+      const bool a_is_finite = (!(a_re_isnan || a_re_isinf || a_im_isnan || a_im_isinf));
+
       complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE> result;
 
-      if(x > BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(0))
+      if(args_is_nan)
+      {
+        result =
+          complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
+          (
+            std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::quiet_NaN(),
+            std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::quiet_NaN()
+          );
+      }
+      else if(x_isinf)
+      {
+        if(a_is_finite)
+        {
+          result =
+            complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
+            (
+              std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::infinity(),
+              std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::infinity()
+            );
+        }
+        else
+        {
+          result =
+            complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
+            (
+              std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::quiet_NaN(),
+              std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::quiet_NaN()
+            );
+        }
+      }
+      else if(x > 0)
       {
         result = std::exp(a * std::log(x));
       }
-      else if(x < BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(0))
+      else if(x < 0)
       {
         using std::acos;
         using std::log;
@@ -531,12 +579,24 @@
       }
       else
       {
-        result =
-          complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
-          (
-            -std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::infinity(),
-            BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(0)
-          );
+        if(a_is_finite)
+        {
+          result =
+            complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
+            (
+              BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(0),
+              BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE(0)
+            );
+        }
+        else
+        {
+          result =
+            complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>
+            (
+              std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::quiet_NaN(),
+              std::numeric_limits<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>::quiet_NaN()
+            );
+        }
       }
 
       return result;
diff --git a/test/float128/test_std_lib.cpp b/test/float128/test_std_lib.cpp
index 661f72c50..d56c7e135 100644
--- a/test/float128/test_std_lib.cpp
+++ b/test/float128/test_std_lib.cpp
@@ -228,6 +228,10 @@ BOOST_AUTO_TEST_CASE( test_main )
    BOOST_CHECK_CLOSE_FRACTION(real(pow(BOOST_FLOATMAX_C(-2.5), cm)), BOOST_FLOATMAX_C(0.0000108384213983921628818716620216475625862428265426558408687988615), 10*tol);
    BOOST_CHECK_CLOSE_FRACTION(imag(pow(BOOST_FLOATMAX_C(-2.5), cm)), BOOST_FLOATMAX_C(-0.0001654255694465738439289663540771804961844719634831318175904326223), 10*tol);
 
+   // Check x^a, where x is zero and a is finite.
+   BOOST_CHECK_EQUAL(real(pow(BOOST_FLOATMAX_C(0.0), cm)), 0);
+   BOOST_CHECK_EQUAL(imag(pow(BOOST_FLOATMAX_C(0.0), cm)), 0);
+
    BOOST_CHECK_CLOSE_FRACTION(real(pow(BOOST_FLOATMAX_C(23.125), cm)), BOOST_FLOATMAX_C(-6.10574617260000071495777483951769228578270070743952693687), 500*tol);
    BOOST_CHECK_CLOSE_FRACTION(imag(pow(BOOST_FLOATMAX_C(23.125), cm)), BOOST_FLOATMAX_C(-2571.59829653692515304089117319850284971907684832627401081405), tol);