Fix complex asin, acos, tan, tanh and log1p mismatches with std

include/xsimd/arch/common/xsimd_common_math.hpp

@@ -1875,11 +1875,13 @@ namespace xsimd
             using real_batch = typename batch_type::real_batch;
             batch_type u = 1 + self;
             batch_type logu = log(u);
-            return select(u == batch_type(1.),
-                          self,
-                          select(u.real() <= real_batch(0.),
-                                 logu,
-                                 logu * self / (u - batch_type(1.))));
+            batch_type res = select(u == batch_type(1.),
+                                    self,
+                                    select(u.real() <= real_batch(0.),
+                                           logu,
+                                           logu * self / (u - batch_type(1.))));
+            // Im(log(1 + z)) = arg(1 + z) has the sign of Im(z). Restore it.
+            return { res.real(), copysign(res.imag(), self.imag()) };
         }
 
         // mod

include/xsimd/arch/common/xsimd_common_trigo.hpp

@@ -46,13 +46,72 @@ namespace xsimd
             x = select(self < batch_type(-0.5), constants::pi<batch_type>() - x, x);
             return select(x_larger_05, x, constants::pio2<batch_type>() - x);
         }
+        namespace detail
+        {
+            // For z = x + iy, returns Re sqrt(1 - z^2) as the real part and |Im asin(z)| as the imaginary part.
+            // With p = sqrt(1 + z) and q = sqrt(1 - z):
+            //   Re sqrt(1 - z^2) = Re(q * p),  |Im asin(z)| = asinh(Im(conj(q) p))
+            // from W. Kahan, "Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit",
+            // in The State of the Art in Numerical Analysis, A. Iserles and M. J. D. Powell (eds.), 1987, pp. 165-211.
+            template <class A, class T>
+            XSIMD_INLINE batch<std::complex<T>, A> asin_kernel(const batch<std::complex<T>, A>& z) noexcept
+            {
+                using batch_type = batch<std::complex<T>, A>;
+                using real_batch = typename batch_type::real_batch;
+                const real_batch one(1.);
+                const real_batch half(0.5);
+                const real_batch x = abs(z.real());
+                const real_batch y = abs(z.imag());
+                const real_batch hy = half * y; // y / 2
+                const real_batch xm1 = x - one; // x - 1
+                const real_batch xp1 = x + one; // x + 1
+
+                // p = sqrt(1 + z) = a + bi, and q = sqrt(1 - z) with a bigger part c and a smaller part d:
+                //   a = sqrt(((1 + x) + sqrt((1 + x)^2 + y^2)) / 2),  b = y / (2a)
+                //   c = sqrt((sqrt((1 - x)^2 + y^2) + |1 - x|) / 2),  d = y / (2c)
+                // 1 + x >= 1, so nothing cancels in a and a >= 1. Using |1 - x| instead of 1 - x adds two
+                // non-negative terms, so nothing cancels in c. For x <= 1, q = c - di (c is the real part); for x > 1,
+                // c and d swap places: q = d - ci.
+                // s = sqrt((1 - x)^2 + y^2). At x = 1 it is y, taken directly: hypot squares y, and y^2 underflows
+                // for tiny y.
+                const real_batch s = select(xm1 == real_batch(0.), y, hypot(xm1, y));
+                const real_batch a = sqrt(half * (hypot(xp1, y) + xp1));
+                const real_batch b = hy / a;
+                const real_batch c = sqrt(half * (s + abs(xm1)));
+                const real_batch d = hy / max(c, constants::smallestposval<real_batch>()); // c is 0 at z = +-1
+                const auto inside = x <= one;
+
+                // Re sqrt(1 - z^2) = Re(q * p).
+                //   x <= 1: (c - di)(a + bi) has real part ac + bd.  x > 1: (d - ci)(a + bi) has real part ad + bc.
+                real_batch re = select(inside, fma(a, c, b * d), fma(a, d, b * c));
+
+                // |Im asin(z)| = asinh(Im(conj(q) p)).
+                //   x <= 1: (c + di)(a + bi) has imaginary part ad + bc.  x > 1: (d + ci)(a + bi) has imaginary part ac + bd.
+                const real_batch sh = select(inside, fma(a, d, b * c), fma(a, c, b * d));
+                real_batch v = asinh(sh);
+
+                // For huge x or y, (1 +- x)^2 + y^2 overflows. There Re sqrt(1 - z^2) = y and
+                // |Im asin(z)| = log(2 sqrt(x^2 + y^2)) up to rounding, as in T. E. Hull, T. F. Fairgrieve and
+                // P. T. P. Tang, "Implementing the Complex Arcsine and Arccosine Functions Using Exception Handling",
+                // ACM Trans. Math. Software 23(3), 1997, pp. 299-335.
+                const auto y_larger = x < y;
+                const real_batch m = select(y_larger, y, x); // max(|x|, |y|)
+                const auto huge = m >= constants::oneotwoeps<real_batch>();
+                if (any(huge))
+                {
+                    const real_batch r = select(y_larger, x / y, y / x); // min(|x|, |y|) / max(|x|, |y|)
+                    re = select(huge, y, re); // Re sqrt(1 - z^2) = |y|
+                    v = select(huge, log(m) + constants::log_2<real_batch>() + half * log1p(r * r), v); // log(2 sqrt(x^2 + y^2)) = log(m) + log(2) + log(1 + r^2) / 2
+                }
+                return { re, v };
+            }
+        }
         template <class A, class T>
         XSIMD_INLINE batch<std::complex<T>, A> acos(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
         {
             using batch_type = batch<std::complex<T>, A>;
-            using real_batch = typename batch_type::real_batch;
-            batch_type tmp = asin(z);
-            return { constants::pio2<real_batch>() - tmp.real(), -tmp.imag() };
+            const batch_type w = detail::asin_kernel(z); // (Re sqrt(1 - z^2), |Im asin(z)|)
+            return { atan2(w.real(), z.real()), copysign(w.imag(), -z.imag()) };
         }
 
         // acosh
@@ -79,9 +138,8 @@ namespace xsimd
         XSIMD_INLINE batch<std::complex<T>, A> acosh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
         {
             using batch_type = batch<std::complex<T>, A>;
-            batch_type w = acos(z);
-            w = batch_type(-w.imag(), w.real());
-            return w;
+            const batch_type w = acos(z);
+            return { abs(w.imag()), copysign(w.real(), z.imag()) };
         }
 
         // asin
@@ -130,20 +188,8 @@ namespace xsimd
         XSIMD_INLINE batch<std::complex<T>, A> asin(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
         {
             using batch_type = batch<std::complex<T>, A>;
-            using real_batch = typename batch_type::real_batch;
-            real_batch x = z.real();
-            real_batch y = z.imag();
-
-            batch_type ct(-y, x);
-            batch_type zz(real_batch(1.) - (x - y) * (x + y), -2 * x * y);
-            zz = log(ct + sqrt(zz));
-            batch_type resg(zz.imag(), -zz.real());
-
-            return select(y == real_batch(0.),
-                          select(fabs(x) > real_batch(1.),
-                                 batch_type(constants::pio2<real_batch>(), real_batch(0.)),
-                                 batch_type(asin(x), real_batch(0.))),
-                          resg);
+            const batch_type w = detail::asin_kernel(z);
+            return { atan2(z.real(), w.real()), copysign(w.imag(), z.imag()) };
         }
 
         // asinh
@@ -354,15 +400,14 @@ namespace xsimd
         XSIMD_INLINE batch<T, A> atan2(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
         {
             using batch_type = batch<T, A>;
+            const batch_type q = select(self == batch_type(0.), batch_type(0.), abs(self / other)); // atan2(+-0, +-0): 0 / 0 would be NaN
 #ifdef __FAST_MATH__
-            const batch_type q = abs(self / other);
             const batch_type q_p = abs(other / self);
 #else
-            const batch_type q = abs(self / other);
             const batch_type q_p = 1. / q;
 #endif
             const batch_type z = detail::kernel_atan(q, q_p);
-            return select(other > batch_type(0.), z, constants::pi<batch_type>() - z) * signnz(self);
+            return select(signnz(other) > batch_type(0.), z, constants::pi<batch_type>() - z) * signnz(self);
         }
 
         // cos
@@ -858,21 +903,13 @@ namespace xsimd
             const batch_type y = detail::tan_eval(xr, test);
             return y ^ bitofsign(self);
         }
+
         template <class A, class T>
         XSIMD_INLINE batch<std::complex<T>, A> tan(batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
         {
             using batch_type = batch<std::complex<T>, A>;
-            using real_batch = typename batch_type::real_batch;
-            real_batch d = cos(2 * z.real()) + cosh(2 * z.imag());
-            real_batch wreal = sin(2 * z.real()) / d;
-            real_batch wimag = sinh(2 * z.imag());
-#ifdef __FAST_MATH__
-            return batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d);
-#else
-            batch_type winf(constants::infinity<real_batch>(), constants::infinity<real_batch>());
-            batch_type wres = select(isinf(wimag), batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d));
-            return select(d == real_batch(0.), winf, wres);
-#endif
+            const batch_type w = tanh(batch_type(-z.imag(), z.real()));
+            return { w.imag(), -w.real() };
         }
 
         // tanh
@@ -978,15 +1015,46 @@ namespace xsimd
             batch_type r = fma(batch_type(-2.), one / (one + exp(x + x)), one);
             return select(test, z, r) ^ bts;
         }
+        // With z = x + iy,
+        //   tanh(z) = ((1 + tan(y)^2) sinh(x) cosh(x) + i tan(y)) / (1 + (1 + tan(y)^2) sinh(x)^2)
+        // from W. Kahan, "Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit" (1987).
+        // It is (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y)) with numerator and denominator divided by
+        // 2 cos(y)^2. For |x| >= 11 (float) or 22 (double) it rounds to
+        //   tanh(z) = sign(x) + i tan(y) / ((1 + tan(y)^2) sinh(x)^2)
+        //
+        // In the code, with x = Re z and y = Im z:
+        //   t = tan(y),  beta = 1 + t^2,  s = sinh(|x|),  c = cosh(|x|) = sqrt(1 + s^2)
+        //   tanh(z) = sign(x) beta s c / (1 + beta s^2) + i t / (1 + beta s^2)
+        //   tanh(z) = sign(x) + i t / (beta s^2)   for |x| >= 11 (float) or 22 (double)
         template <class A, class T>
         XSIMD_INLINE batch<std::complex<T>, A> tanh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
         {
-            using real_batch = typename batch<std::complex<T>, A>::real_batch;
-            auto x = z.real();
-            auto y = z.imag();
-            real_batch two(2);
-            auto d = cosh(two * x) + cos(two * y);
-            return { sinh(two * x) / d, sin(two * y) / d };
+            using batch_type = batch<std::complex<T>, A>;
+            using real_batch = typename batch_type::real_batch;
+            const real_batch one(1.);
+            const real_batch x = z.real();
+            const real_batch ax = abs(x);
+
+            // t and beta = 1 / cos(y)^2 >= 1.
+            const real_batch t = tan(z.imag());
+            const real_batch beta = fma(t, t, one);
+
+            // s and c. |x| is capped so that s stays finite: with -ffast-math, sinh of a large |x| is
+            // garbage instead of inf, and 1 / inf is NaN.
+            const real_batch cap(std::is_same_v<T, float> ? T(64) : T(512));
+            const real_batch s = sinh(select(ax > cap, cap, ax));
+            const real_batch c = sqrt(fma(s, s, one));
+
+            // Main formula. All terms are >= 0 and the denominator is >= 1: no cancellation, no division by 0.
+            const real_batch inv = one / fma(beta * s, s, one);
+            const real_batch re = beta * s * c * inv;
+            const real_batch im = t * inv;
+
+            // Large |x|. There 1 + beta s^2 rounds to beta s^2 and c / s rounds to 1, so re = sign(x) and
+            // im = t / (beta s^2), computed as (t / beta) (1 / s) (1 / s) because s^2 can overflow.
+            const auto large = ax >= real_batch(std::is_same_v<T, float> ? T(11) : T(22));
+            const real_batch inv_s = one / s;
+            return { copysign(select(large, one, re), x), select(large, (t / beta) * inv_s * inv_s, im) };
         }
 
     }

include/xsimd/arch/xsimd_scalar.hpp

@@ -701,7 +701,9 @@ namespace xsimd
         {
             using T = typename C::value_type;
             C u = C(1.) + val;
-            return u == C(1.) ? val : (u.real() <= T(0.) ? log(u) : log(u) * val / (u - C(1.)));
+            C res = u == C(1.) ? val : (u.real() <= T(0.) ? log(u) : log(u) * val / (u - C(1.)));
+            // Im(log(1 + z)) has the sign of Im(z), which the operations above lose for an imaginary -0.
+            return C(res.real(), std::copysign(res.imag(), val.imag()));
         }
     }
 

test/test_complex_exponential.cpp

@@ -213,4 +213,27 @@ TEST_CASE_TEMPLATE("[complex exponential]", B, BATCH_COMPLEX_TYPES)
         Test.test_sign();
     }
 }
+
+// Cases of issue 1432 that the fix changed, as in test_complex_trigonometric.cpp.
+TEST_CASE_TEMPLATE("[complex exponential issue 1432]", B, BATCH_COMPLEX_TYPES)
+{
+    using xsimd::constants::log_2;
+    using xsimd::constants::pi;
+
+    std::vector<std::complex<double>> inputs = {
+        { 1, -0. }, // imaginary part -0
+        { -2, -0. }, // the -0 sets the sign of Im log1p(z)
+    };
+    std::vector<std::complex<double>> expected = {
+        { log_2<double>(), -0. },
+        { 0., -pi<double>() },
+    };
+    std::vector<double> cond = { 1, 1 };
+#ifndef __FAST_MATH__ // no infinities with fast-math
+    inputs.push_back({ -1, 0. });
+    expected.push_back({ -std::numeric_limits<double>::infinity(), 0. });
+    cond.push_back(1);
+#endif
+    CHECK_COMPLEX_CASES(B, log1p, inputs, expected, cond);
+}
 #endif

test/test_complex_trigonometric.cpp

@@ -229,4 +229,116 @@ TEST_CASE_TEMPLATE("[complex trigonometric]", B, BATCH_COMPLEX_TYPES)
         Test.test_atan();
     }
 }
+
+// Cases of issue 1432 that the fix changed: each gave a wrong value or a wrong sign of zero before it. The expected
+// values are the exact f(z) rounded to double; cond is max(1, |z f'(z) / f(z)|) (see detail::check_complex_cases).
+TEST_CASE_TEMPLATE("[complex trigonometric issue 1432]", B, BATCH_COMPLEX_TYPES)
+{
+    using xsimd::constants::pi;
+    using xsimd::constants::pio2;
+
+    SUBCASE("asin")
+    {
+        const std::vector<std::complex<double>> inputs = {
+            { -2, -0. }, // the -0 sets the sign of Im asin(z)
+            { -0., -0. }, // signed zeros
+            { 1, 800 }, // large imaginary part
+            { 10000, 1 }, // large real part
+        };
+        const std::vector<std::complex<double>> expected = {
+            { -pio2<double>(), -1.3169578969248168 },
+            { -0., -0. },
+            { 0.0012499983724006398, 7.377760080100202 },
+            { 1.57069632679473, 9.903487555036127 },
+        };
+        const std::vector<double> cond = { 1, 1, 1, 1 };
+        CHECK_COMPLEX_CASES(B, asin, inputs, expected, cond);
+    }
+
+    SUBCASE("acos")
+    {
+        const std::vector<std::complex<double>> inputs = {
+            { -2, -0. }, // the -0 sets the sign of Im acos(z)
+            { 0., -0. }, // signed zeros
+            { 1, 800 }, // large imaginary part
+            { 10000, 1 }, // large real part
+        };
+        const std::vector<std::complex<double>> expected = {
+            { pi<double>(), 1.3169578969248168 },
+            { pio2<double>(), 0. },
+            { 1.569546328422496, -7.377760080100202 },
+            { 0.00010000000016666667, -9.903487555036127 },
+        };
+        const std::vector<double> cond = { 1, 1, 1, 1 };
+        CHECK_COMPLEX_CASES(B, acos, inputs, expected, cond);
+    }
+
+    SUBCASE("tan")
+    {
+        const std::vector<std::complex<double>> inputs = {
+            { 1.5707f, 0. }, // near the pole pi/2 = 1.570796
+            { 2, -500 }, // large imaginary part
+        };
+        const std::vector<std::complex<double>> expected = {
+            { 10386.654119664716, 0. },
+            { -0., -1. },
+        };
+        const std::vector<double> cond = { 1.63e4, 1 };
+        CHECK_COMPLEX_CASES(B, tan, inputs, expected, cond);
+    }
+
+    SUBCASE("asinh")
+    {
+        const std::vector<std::complex<double>> inputs = {
+            { -0., -3 }, // the -0 sets the sign of Re asinh(z)
+            { -0., -0. }, // signed zeros
+            { 2, 0. }, // real axis
+            { 1, -500 }, // large imaginary part
+            { 10000, 1 }, // large real part
+        };
+        const std::vector<std::complex<double>> expected = {
+            { -1.762747174039086, -pio2<double>() },
+            { -0., -0. },
+            { 1.4436354751788103, 0. },
+            { 6.907756278988637, -1.5687963254615769 },
+            { 9.903487560036128, 9.999999916666668e-05 },
+        };
+        const std::vector<double> cond = { 1, 1, 1, 1, 1 };
+        CHECK_COMPLEX_CASES(B, asinh, inputs, expected, cond);
+    }
+
+    SUBCASE("acosh")
+    {
+        const std::vector<std::complex<double>> inputs = {
+            { -4, -3 }, // x < 0, y < 0
+            { -2, 0. }, // the +0 sets the sign of Im acosh(z)
+            { 1, -0. }, // signed zeros
+            { 1, 800 }, // large imaginary part
+            { 10000, 1 }, // large real part
+        };
+        const std::vector<std::complex<double>> expected = {
+            { 2.2999140408792695, -2.488413180146375 },
+            { 1.3169578969248168, pi<double>() },
+            { 0., -0. },
+            { 7.377760080100202, 1.569546328422496 },
+            { 9.903487555036127, 0.00010000000016666667 },
+        };
+        const std::vector<double> cond = { 1, 1, 1, 1, 1 };
+        CHECK_COMPLEX_CASES(B, acosh, inputs, expected, cond);
+    }
+
+    SUBCASE("tanh")
+    {
+        const std::vector<std::complex<double>> inputs = {
+            { -400, 2 }, // large real part
+            { 0., -10.9958f }, // near the pole -7 pi/2 i = -10.995574 i
+        };
+        const std::vector<std::complex<double>> expected = {
+            { -1., -0. },
+            { 0., 4430.056602996832 },
+        };
+        const std::vector<double> cond = { 1, 4.87e4 };
+        CHECK_COMPLEX_CASES(B, tanh, inputs, expected, cond);
+    }
+}
 #endif