[Bug libstdc++/107538] std::pow(10, std::complex(NaN, 1)) aborts with -D_GLIBCXX_ASSERTIONS

[redi at gcc dot gnu.org via Gcc-bugs]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107538

--- Comment #3 from Jonathan Wakely <redi at gcc dot gnu.org> ---
(In reply to Jonathan Wakely from comment #0)
> #define _GLIBCXX_ASSERTIONS
> #include <complex>
> #include <limits>
> 
> int main()
> {
>   double nan = std::numeric_limits<double>::quiet_NaN();
>   std::pow(10, std::complex<double>(nan, 1));

Actually maybe the assertion is OK. This violates [complex.numbers] p3:

"If the result of a function is not mathematically defined or not in the range
of representable values for its type, the behavior is undefined."

10^(nan + 1i) is not mathematically defined, an assertion seems OK.

If we want to prevent the assertion here we could do:

--- a/libstdc++-v3/include/std/complex
+++ b/libstdc++-v3/include/std/complex
@@ -699,12 +699,19 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
        && !_GLIBCXX_FAST_MATH>::_S_do_it(__z);
     }

+  template<typename _Tp>
+    inline complex<_Tp>
+    __polar(const _Tp& __rho, const _Tp& __theta)
+    { return complex<_Tp>(__rho * cos(__theta), __rho * sin(__theta)); }
+
   template<typename _Tp>
     inline complex<_Tp>
     polar(const _Tp& __rho, const _Tp& __theta)
     {
+      // _GLIBCXX_RESOLVE_LIB_DEFECTS
+      // 2459. std::polar should require a non-negative rho
       __glibcxx_assert( __rho >= 0 );
-      return complex<_Tp>(__rho * cos(__theta), __rho * sin(__theta));
+      return std::__polar<_Tp>(__rho, __theta);
     }

   template<typename _Tp>
@@ -778,7 +785,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
   template<typename _Tp>
     inline complex<_Tp>
     __complex_exp(const complex<_Tp>& __z)
-    { return std::polar<_Tp>(exp(__z.real()), __z.imag()); }
+    { return std::__polar<_Tp>(exp(__z.real()), __z.imag()); }

 #if _GLIBCXX_USE_C99_COMPLEX
   inline __complex__ float
@@ -1038,7 +1045,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
         return pow(__x.real(), __y);

       complex<_Tp> __t = std::log(__x);
-      return std::polar<_Tp>(exp(__y * __t.real()), __y * __t.imag());
+      return std::__polar<_Tp>(exp(__y * __t.real()), __y * __t.imag());
     }

   template<typename _Tp>
@@ -1075,8 +1082,8 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
     inline complex<_Tp>
     pow(const _Tp& __x, const complex<_Tp>& __y)
     {
-      return __x > _Tp() ? std::polar<_Tp>(pow(__x, __y.real()),
-                                          __y.imag() * log(__x))
+      return __x > _Tp() ? std::__polar<_Tp>(pow(__x, __y.real()),
+                                            __y.imag() * log(__x))
                         : std::pow(complex<_Tp>(__x), __y);
     }