[libcxx] Improve accuracy of complex exp

[nurmukhametov]

When computing exp(x + iy) with large x and small y, the naive computation exp(x) * (cos(y) + i*sin(y)) produces inf when exp(x) overflows, even though the result is finite when y is sufficiently small. Fix this with exp(x/2) * (cos(y) + i*sin(y)) * exp(x/2)

[llvmbot]

@llvm/pr-subscribers-libcxx

Author: Aleksei Nurmukhametov (nurmukhametov)

Changes

When computing exp(x + iy) with large x and small y, the naive computation exp(x) * (cos(y) + i*sin(y)) produces inf when exp(x) overflows, even though the result is finite when y is sufficiently small. Fix this with exp(x/2) * (cos(y) + i*sin(y)) * exp(x/2)

---

Full diff: https://github.com/llvm/llvm-project/pull/165254.diff

2 Files Affected:

• (modified) libcxx/include/complex (+4)
• (added) libcxx/test/libcxx/numerics/complex.number/exp.pass.cpp (+44)

diff --git a/libcxx/include/complex b/libcxx/include/complex
index d8ec3d95c10ed..c749156e5d4f5 100644
--- a/libcxx/include/complex
+++ b/libcxx/include/complex
@@ -1091,6 +1091,10 @@ _LIBCPP_HIDE_FROM_ABI complex<_Tp> exp(const complex<_Tp>& __x) {
     }
   }
   _Tp __e = std::exp(__x.real());
+  if (std::isinf(__e)) {
+    _Tp __e2 = std::exp(__x.real() * _Tp(0.5));
+    return complex<_Tp>(__e2 * std::cos(__i) * __e2, __e2 * std::sin(__i) * __e2);
+  }
   return complex<_Tp>(__e * std::cos(__i), __e * std::sin(__i));
 }
 
diff --git a/libcxx/test/libcxx/numerics/complex.number/exp.pass.cpp b/libcxx/test/libcxx/numerics/complex.number/exp.pass.cpp
new file mode 100644
index 0000000000000..74c4a1887a8ec
--- /dev/null
+++ b/libcxx/test/libcxx/numerics/complex.number/exp.pass.cpp
@@ -0,0 +1,44 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+// <complex>
+
+// template<class T>
+//   complex<T>
+//   exp(const complex<T>& x);
+//
+// Tests for libc++-specific overflow handling behavior in complex exponential.
+// These tests validate implementation-specific handling of edge cases where
+// exp(real_part) overflows but the result should still be well-defined.
+
+#include <complex>
+#include <cassert>
+#include <cmath>
+
+#include "test_macros.h"
+
+template <class T>
+void test_overflow_case() {
+    typedef std::complex<T> C;
+
+    // In this case, the overflow of exp(real_part) is compensated when
+    // sin(imag_part) is close to zero, resulting in a finite imaginary part.
+    C z(T(90.0238094), T(5.900613e-39));
+    C result = std::exp(z);
+
+    assert(std::isinf(result.real()));
+    assert(result.real() > 0);
+
+    assert(std::isfinite(result.imag()));
+    assert(std::abs(result.imag() - T(7.3746)) < T(1.0));
+}
+
+int main(int, char**) {
+    test_overflow_case<float>();
+    return 0;
+}

[github-actions]

✅ With the latest revision this PR passed the C/C++ code formatter.

[philnik777]

I think an approach similar to #99677 would be better.

[nurmukhametov]

I think an approach similar to #99677 would be better.

I didn't know about it. I can do the same, but I wonder why it hasn't been merged yet.

[philnik777]

I think an approach similar to #99677 would be better.

I didn't know about it. I can do the same, but I wonder why it hasn't been merged yet.

IIRC it's mostly due to me having no idea what I'm doing with floating point numbers and the CI complaining in some place. Anybody who know anything about them probably knows what to do. (So if you want to after this patch, feel free to take a look that the rest of the functions. That would be very much appreciated.)

[nurmukhametov]

I think an approach similar to #99677 would be better.

I have created this PR #165457 with this approach. However, I don't really understand why you preserved partially the initial implementation in #99677.

[nurmukhametov]

I think an approach similar to #99677 would be better.

This approach (calling implementations of underlying libm) doesn't fix the initial issue (the accuracy of complex exponential implementation) the current PR addressing. You can see it in #165457 for example for CI jobs with picolibc.

[philnik777]

I think an approach similar to #99677 would be better.

This approach (calling implementations of underlying libm) doesn't fix the initial issue (the accuracy of complex exponential implementation) the current PR addressing. You can see it in #165457 for example for CI jobs with picolibc.

I don't think that's a huge problem. We're really not in the business of implementing high quality math functions (see <complex>). That's the kind of stuff we should leave to the libc. We wouldn't implement any other math functions just because some libc does a sub-par job. If your libc doesn't provide a high quality enough implementation for your purposes you can file a bug against that or use a LLVM libc shim instead (once they implement the complex family of functions).

[nurmukhametov]

I think an approach similar to #99677 would be better.

This approach (calling implementations of underlying libm) doesn't fix the initial issue (the accuracy of complex exponential implementation) the current PR addressing. You can see it in #165457 for example for CI jobs with picolibc.

I don't think that's a huge problem. We're really not in the business of implementing high quality math functions (see <complex>). That's the kind of stuff we should leave to the libc. We wouldn't implement any other math functions just because some libc does a sub-par job. If your libc doesn't provide a high quality enough implementation for your purposes you can file a bug against that or use a LLVM libc shim instead (once they implement the complex family of functions).

I understand this point, but it is not so currently. At the moment, libcxx has its own implementations.