<complex>: Fix exp and polar on NaN/infinity input (#1584)

Co-authored-by: Stephan T. Lavavej <stl@microsoft.com>

stl/inc/complex

@@ -1813,12 +1813,46 @@ _NODISCARD complex<_Ty> cosh(const complex<_Ty>& _Left) {
         _Ctraits<_Ty>::_Sinh(real(_Left), _Ctraits<_Ty>::sin(imag(_Left))));
 }
 
+template <class _Ty>
+_NODISCARD complex<_Ty> _Polar_positive_nan_inf_zero_rho(const _Ty& _Rho, const _Ty& _Theta) { // _Rho is +NaN/+Inf/+0
+    if (_Ctraits<_Ty>::_Isnan(_Theta) || _Ctraits<_Ty>::_Isinf(_Theta)) { // _Theta is NaN/Inf
+        if (_Ctraits<_Ty>::_Isinf(_Rho)) {
+            return complex<_Ty>(_Rho, _Ctraits<_Ty>::sin(_Theta)); // (Inf, NaN/Inf)
+        } else {
+            return complex<_Ty>(_Rho, _Ctraits<_Ty>::_Copysign(_Rho, _Theta)); // (NaN/0, NaN/Inf)
+        }
+    } else if (_Theta == _Ty{0}) { // _Theta is zero
+        return complex<_Ty>(_Rho, _Theta); // (NaN/Inf/0, 0)
+    } else { // _Theta is finite non-zero
+        // (NaN/Inf/0, finite non-zero)
+        return complex<_Ty>(_Rho * _Ctraits<_Ty>::cos(_Theta), _Rho * _Ctraits<_Ty>::sin(_Theta));
+    }
+}
+
 template <class _Ty>
 _NODISCARD complex<_Ty> exp(const complex<_Ty>& _Left) {
-    _Ty _Real(real(_Left)), _Imag(real(_Left));
-    _Ctraits<_Ty>::_Exp(&_Real, _Ctraits<_Ty>::cos(imag(_Left)), 0);
-    _Ctraits<_Ty>::_Exp(&_Imag, _Ctraits<_Ty>::sin(imag(_Left)), 0);
-    return complex<_Ty>(_Real, _Imag);
+    const _Ty _Log_rho = _STD real(_Left);
+    const _Ty _Theta   = _STD imag(_Left);
+
+    if (!_Ctraits<_Ty>::_Isnan(_Log_rho) && !_Ctraits<_Ty>::_Isinf(_Log_rho)) { // real component is finite
+        _Ty _Real = _Log_rho;
+        _Ty _Imag = _Log_rho;
+        _Ctraits<_Ty>::_Exp(&_Real, _Ctraits<_Ty>::cos(_Theta), 0);
+        _Ctraits<_Ty>::_Exp(&_Imag, _Ctraits<_Ty>::sin(_Theta), 0);
+        return complex<_Ty>(_Real, _Imag);
+    }
+
+    // real component is NaN/Inf
+    // return polar(exp(re), im)
+    if (_Ctraits<_Ty>::_Isinf(_Log_rho)) {
+        if (_Log_rho < _Ty{0}) {
+            return _Polar_positive_nan_inf_zero_rho(_Ty{0}, _Theta); // exp(-Inf) = +0
+        } else {
+            return _Polar_positive_nan_inf_zero_rho(_Log_rho, _Theta); // exp(+Inf) = +Inf
+        }
+    } else {
+        return _Polar_positive_nan_inf_zero_rho(_Ctraits<_Ty>::_Abs(_Log_rho), _Theta); // exp(NaN) = +NaN
+    }
 }
 
 template <class _Ty>
@@ -2057,7 +2091,16 @@ _NODISCARD _CONSTEXPR20 _Ty norm(const complex<_Ty>& _Left) { // return squared
 
 template <class _Ty>
 _NODISCARD complex<_Ty> polar(const _Ty& _Rho, const _Ty& _Theta) { // return _Rho * exp(i * _Theta) as complex
-    return complex<_Ty>(_Rho * _Ctraits<_Ty>::cos(_Theta), _Rho * _Ctraits<_Ty>::sin(_Theta));
+    if (!_Ctraits<_Ty>::_Isnan(_Rho) && !_Ctraits<_Ty>::_Isinf(_Rho) && _Rho != _Ty{0}) { // _Rho is finite non-zero
+        return complex<_Ty>(_Rho * _Ctraits<_Ty>::cos(_Theta), _Rho * _Ctraits<_Ty>::sin(_Theta));
+    }
+
+    // _Rho is NaN/Inf/0
+    if (_Ctraits<_Ty>::_Signbit(_Rho)) {
+        return -_Polar_positive_nan_inf_zero_rho(-_Rho, _Theta);
+    } else {
+        return _Polar_positive_nan_inf_zero_rho(_Rho, _Theta);
+    }
 }
 
 template <class _Ty>

stl/src/xexp.cpp

@@ -17,11 +17,14 @@ static const double invln2 = 1.4426950408889634073599246810018921;
 
 _CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _Exp(
     double* px, double y, short eoff) { // compute y * e^(*px), (*px) finite, |y| not huge
-    if (*px < -hugexp || y == 0.0) { // certain underflow
-        *px = 0.0;
+    if (y == 0.0) { // zero
+        *px = y;
+        return 0;
+    } else if (*px < -hugexp) { // certain underflow
+        *px = _Xfe_underflow(y);
         return 0;
     } else if (hugexp < *px) { // certain overflow
-        *px = _Inf._Double;
+        *px = _Xfe_overflow(y);
         return _INFCODE;
     } else { // xexp won't overflow
         double g   = *px * invln2;
@@ -39,7 +42,21 @@ _CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _Exp(
             *px = (w + g) / (w - g) * 2.0 * y;
             --xexp;
         }
-        return _Dscale(px, static_cast<long>(xexp) + eoff);
+
+        const short result_code = _Dscale(px, static_cast<long>(xexp) + eoff);
+
+        switch (result_code) {
+        case 0:
+            *px = _Xfe_underflow(y);
+            break;
+        case _INFCODE:
+            *px = _Xfe_overflow(y);
+            break;
+        default:
+            break;
+        }
+
+        return result_code;
     }
 }
 

stl/src/xfexp.cpp

@@ -16,11 +16,14 @@ static const float invln2 = 1.4426950408889634073599246810018921F;
 
 _CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _FExp(
     float* px, float y, short eoff) { // compute y * e^(*px), (*px) finite, |y| not huge
-    if (*px < -hugexp || y == 0.0F) { // certain underflow
-        *px = 0.0F;
+    if (y == 0.0F) { // zero
+        *px = y;
+        return 0;
+    } else if (*px < -hugexp) { // certain underflow
+        *px = _Xfe_underflow(y);
         return 0;
     } else if (hugexp < *px) { // certain overflow
-        *px = _FInf._Float;
+        *px = _Xfe_overflow(y);
         return _INFCODE;
     } else { // xexp won't overflow
         float g    = *px * invln2;
@@ -38,7 +41,21 @@ _CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _FExp(
             *px = (w + g) / (w - g) * 2.0F * y;
             --xexp;
         }
-        return _FDscale(px, static_cast<long>(xexp) + eoff);
+
+        const short result_code = _FDscale(px, static_cast<long>(xexp) + eoff);
+
+        switch (result_code) {
+        case 0:
+            *px = _Xfe_underflow(y);
+            break;
+        case _INFCODE:
+            *px = _Xfe_overflow(y);
+            break;
+        default:
+            break;
+        }
+
+        return result_code;
     }
 }
 

stl/src/xlexp.cpp

@@ -18,11 +18,14 @@ static const long double invln2 = 1.4426950408889634073599246810018921L;
 
 _CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _LExp(
     long double* px, long double y, short eoff) { // compute y * e^(*px), (*px) finite, |y| not huge
-    if (*px < -hugexp || y == 0.0L) { // certain underflow
-        *px = 0.0L;
+    if (y == 0.0L) { // zero
+        *px = y;
+        return 0;
+    } else if (*px < -hugexp) { // certain underflow
+        *px = _Xfe_underflow(y);
         return 0;
     } else if (hugexp < *px) { // certain overflow
-        *px = _LInf._Long_double;
+        *px = _Xfe_overflow(y);
         return _INFCODE;
     } else { // xexp won't overflow
         long double g = *px * invln2;
@@ -40,7 +43,21 @@ _CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _LExp(
             *px = (w + g) / (w - g) * 2.0L * y;
             --xexp;
         }
-        return _LDscale(px, static_cast<long>(xexp) + eoff);
+
+        const short result_code = _LDscale(px, static_cast<long>(xexp) + eoff);
+
+        switch (result_code) {
+        case 0:
+            *px = _Xfe_underflow(y);
+            break;
+        case _INFCODE:
+            *px = _Xfe_overflow(y);
+            break;
+        default:
+            break;
+        }
+
+        return result_code;
     }
 }
 

stl/src/xmath.hpp

@@ -5,6 +5,7 @@
 #define _XMATH
 #include <cerrno>
 #include <cmath>
+#include <limits>
 #include <ymath.h>
 
 // macros for _Feraise argument
@@ -165,4 +166,27 @@ long double* _LXp_sqrtx(long double*, int, long double*);
 
 _END_EXTERN_C_UNLESS_PURE
 
+// raise IEEE 754 exceptions
+#ifndef _M_CEE_PURE
+#pragma float_control(except, on, push)
+#endif
+
+template <typename T>
+_NODISCARD T _Xfe_overflow(const T sign) noexcept {
+    static_assert(_STD is_floating_point_v<T>, "Expected is_floating_point_v<T>.");
+    constexpr T huge = _STD numeric_limits<T>::max();
+    return _STD copysign(huge, sign) * huge;
+}
+
+template <typename T>
+_NODISCARD T _Xfe_underflow(const T sign) noexcept {
+    static_assert(_STD is_floating_point_v<T>, "Expected is_floating_point_v<T>.");
+    constexpr T tiny = _STD numeric_limits<T>::min();
+    return _STD copysign(tiny, sign) * tiny;
+}
+
+#ifndef _M_CEE_PURE
+#pragma float_control(pop)
+#endif
+
 #endif // _XMATH

tests/libcxx/expected_results.txt

@@ -557,6 +557,11 @@ std/numerics/complex.number/cmplx.over/conj.pass.cpp:0 FAIL
 std/numerics/complex.number/cmplx.over/pow.pass.cpp:0 FAIL
 std/numerics/complex.number/cmplx.over/proj.pass.cpp:0 FAIL
 
+# Assertion failed: c == NaN || c == non_zero_nan
+# Testing input values outside the range of [complex.value.ops]/9
+# libc++ handles those input values differently
+std/numerics/complex.number/complex.value.ops/polar.pass.cpp FAIL
+
 # Assertion failed: std::abs(skew - x_skew) < 0.01
 # Random number generation test with too strict pass criteria (test8 failure probability ~= 0.04)
 std/numerics/rand/rand.dis/rand.dist.samp/rand.dist.samp.pconst/eval.pass.cpp FAIL
@@ -720,7 +725,6 @@ std/numerics/complex.number/complex.transcendentals/asinh.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/atanh.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/cos.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/cosh.pass.cpp FAIL
-std/numerics/complex.number/complex.transcendentals/exp.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/log10.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/pow_complex_complex.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/pow_complex_scalar.pass.cpp FAIL
@@ -729,7 +733,6 @@ std/numerics/complex.number/complex.transcendentals/sin.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/sinh.pass.cpp FAIL
 std/numerics/complex.number/complex.transcendentals/tanh.pass.cpp FAIL
 std/numerics/complex.number/complex.value.ops/norm.pass.cpp FAIL
-std/numerics/complex.number/complex.value.ops/polar.pass.cpp FAIL
 
 # Not yet analyzed, likely STL bugs. Many various assertions.
 std/localization/locale.categories/category.ctype/facet.ctype.special/facet.ctype.char.statics/classic_table.pass.cpp FAIL

tests/libcxx/skipped_tests.txt

@@ -557,6 +557,11 @@ numerics\complex.number\cmplx.over\conj.pass.cpp
 numerics\complex.number\cmplx.over\pow.pass.cpp
 numerics\complex.number\cmplx.over\proj.pass.cpp
 
+# Assertion failed: c == NaN || c == non_zero_nan
+# Testing input values outside the range of [complex.value.ops]/9
+# libc++ handles those input values differently
+numerics\complex.number\complex.value.ops\polar.pass.cpp
+
 # Assertion failed: std::abs(skew - x_skew) < 0.01
 # Random number generation test with too strict pass criteria (test8 failure probability ~= 0.04)
 numerics\rand\rand.dis\rand.dist.samp\rand.dist.samp.pconst\eval.pass.cpp
@@ -720,7 +725,6 @@ numerics\complex.number\complex.transcendentals\asinh.pass.cpp
 numerics\complex.number\complex.transcendentals\atanh.pass.cpp
 numerics\complex.number\complex.transcendentals\cos.pass.cpp
 numerics\complex.number\complex.transcendentals\cosh.pass.cpp
-numerics\complex.number\complex.transcendentals\exp.pass.cpp
 numerics\complex.number\complex.transcendentals\log10.pass.cpp
 numerics\complex.number\complex.transcendentals\pow_complex_complex.pass.cpp
 numerics\complex.number\complex.transcendentals\pow_complex_scalar.pass.cpp
@@ -729,7 +733,6 @@ numerics\complex.number\complex.transcendentals\sin.pass.cpp
 numerics\complex.number\complex.transcendentals\sinh.pass.cpp
 numerics\complex.number\complex.transcendentals\tanh.pass.cpp
 numerics\complex.number\complex.value.ops\norm.pass.cpp
-numerics\complex.number\complex.value.ops\polar.pass.cpp
 
 # Not yet analyzed, likely STL bugs. Many various assertions.
 localization\locale.categories\category.ctype\facet.ctype.special\facet.ctype.char.statics\classic_table.pass.cpp