[SYCL][libdevice] improve complex trig and pi ulp

libdevice/device_complex.h

@@ -10,6 +10,18 @@
 
 #include "device.h"
 
+#if __cplusplus >= 202002L && __has_include(<version>)
+#include <version> // defines __cpp_lib_math_constants
+#endif
+// Correctly rounded pi for the inverse trig/hyperbolic functions; the
+// atan2(+0, -0) idiom from libc++ is 1 ulp low for float on Intel GPUs.
+#if __cpp_lib_math_constants
+#include <numbers>
+#define __DEVICELIB_PI std::numbers::pi_v<double>
+#else
+#define __DEVICELIB_PI 3.14159265358979323846
+#endif
+
 #if defined(__SPIR__) || defined(__SPIRV__)
 
 // TODO: This needs to be more robust.

libdevice/fallback-complex-fp64.hpp

@@ -191,7 +191,7 @@ static inline double __complex__ __sqr(double __complex__ z) {
 static inline double __complex__ __devicelib_cacos(double __complex__ z) {
   double z_real = __devicelib_creal(z);
   double z_imag = __devicelib_cimag(z);
-  const double __pi(__spirv_ocl_atan2(+0.0, -0.0));
+  const double __pi(__DEVICELIB_PI);
   if (__spirv_IsInf(z_real)) {
     if (__spirv_IsNan(z_imag))
       return CMPLX(z_imag, z_real);
@@ -213,19 +213,25 @@ static inline double __complex__ __devicelib_cacos(double __complex__ z) {
     return CMPLX(__pi / 2.0, -z_imag);
   if (z_real == 0 && (z_imag == 0 || __spirv_IsNan(z_imag)))
     return CMPLX(__pi / 2.0, -z_imag);
+  // acos(z) = pi - acos(-z) for real(z) < 0; evaluating on the right
+  // half-plane keeps z + sqrt(z^2-1) from cancelling.  Im(acos) keeps the
+  // sign opposite to imag(z) either way.
+  const bool refl = __spirv_SignBitSet(z_real);
+  double __complex__ zr = refl ? -z : z;
   double __complex__ w =
-      __devicelib_clog(z + __devicelib_csqrt(__sqr(z) - 1.0));
+      __devicelib_clog(zr + __devicelib_csqrt(__sqr(zr) - 1.0));
+  double re = __spirv_ocl_fabs(__devicelib_cimag(w));
+  if (refl)
+    re = __pi - re;
   if (__spirv_SignBitSet(z_imag))
-    return CMPLX(__spirv_ocl_fabs(__devicelib_cimag(w)),
-                 __spirv_ocl_fabs(__devicelib_creal(w)));
-  return CMPLX(__spirv_ocl_fabs(__devicelib_cimag(w)),
-               -__spirv_ocl_fabs(__devicelib_creal(w)));
+    return CMPLX(re, __spirv_ocl_fabs(__devicelib_creal(w)));
+  return CMPLX(re, -__spirv_ocl_fabs(__devicelib_creal(w)));
 }
 
 static inline double __complex__ __devicelib_casinh(double __complex__ z) {
   double z_real = __devicelib_creal(z);
   double z_imag = __devicelib_cimag(z);
-  const double __pi(__spirv_ocl_atan2(+0.0, -0.0));
+  const double __pi(__DEVICELIB_PI);
   if (__spirv_IsInf(z_real)) {
     if (__spirv_IsNan(z_imag))
       return z;
@@ -243,8 +249,11 @@ static inline double __complex__ __devicelib_casinh(double __complex__ z) {
   if (__spirv_IsInf(z_imag))
     return CMPLX(__spirv_ocl_copysign(z_imag, z_real),
                  __spirv_ocl_copysign(__pi / 2.0, z_imag));
+  // asinh is odd: evaluate on the right half-plane, where sqrt(z^2+1) ~ +z and
+  // the sum cannot cancel, then take the signs from the original argument.
+  double __complex__ zr = __spirv_SignBitSet(z_real) ? -z : z;
   double __complex__ w =
-      __devicelib_clog(z + __devicelib_csqrt(__sqr(z) + 1.0));
+      __devicelib_clog(zr + __devicelib_csqrt(__sqr(zr) + 1.0));
   return CMPLX(__spirv_ocl_copysign(__devicelib_creal(w), z_real),
                __spirv_ocl_copysign(__devicelib_cimag(w), z_imag));
 }
@@ -258,7 +267,7 @@ static inline double __complex__ __devicelib_casin(double __complex__ z) {
 static inline double __complex__ __devicelib_cacosh(double __complex__ z) {
   double z_real = __devicelib_creal(z);
   double z_imag = __devicelib_cimag(z);
-  const double __pi(__spirv_ocl_atan2(+0.0, -0.0));
+  const double __pi(__DEVICELIB_PI);
   if (__spirv_IsInf(z_real)) {
     if (__spirv_IsNan(z_imag))
       return CMPLX(__spirv_ocl_fabs(z_real), z_imag);
@@ -280,16 +289,23 @@ static inline double __complex__ __devicelib_cacosh(double __complex__ z) {
   if (__spirv_IsInf(z_imag))
     return CMPLX(__spirv_ocl_fabs(z_imag),
                  __spirv_ocl_copysign(__pi / 2.0, z_imag));
+  // acosh(z) = acosh(-z) + i*copysign(pi, imag(z)) for real(z) < 0; evaluating
+  // on the right half-plane keeps z + sqrt(z^2-1) from cancelling.
+  const bool refl = __spirv_SignBitSet(z_real);
+  double __complex__ zr = refl ? -z : z;
   double __complex__ w =
-      __devicelib_clog(z + __devicelib_csqrt(__sqr(z) - 1.0));
+      __devicelib_clog(zr + __devicelib_csqrt(__sqr(zr) - 1.0));
+  double im = __spirv_ocl_fabs(__devicelib_cimag(w));
+  if (refl)
+    im = __pi - im;
   return CMPLX(__spirv_ocl_copysign(__devicelib_creal(w), 0.0),
-               __spirv_ocl_copysign(__devicelib_cimag(w), z_imag));
+               __spirv_ocl_copysign(im, z_imag));
 }
 
 static inline double __complex__ __devicelib_catanh(double __complex__ z) {
   double z_real = __devicelib_creal(z);
   double z_imag = __devicelib_cimag(z);
-  const double __pi(__spirv_ocl_atan2(+0.0, -0.0));
+  const double __pi(__DEVICELIB_PI);
   if (__spirv_IsInf(z_imag))
     return CMPLX(__spirv_ocl_copysign(0.0, z_real),
                  __spirv_ocl_copysign(__pi / 2.0, z_imag));

libdevice/fallback-complex.hpp

@@ -184,7 +184,7 @@ static inline float __complex__ __sqrf(float __complex__ z) {
 static inline float __complex__ __devicelib_cacosf(float __complex__ z) {
   float z_real = __devicelib_crealf(z);
   float z_imag = __devicelib_cimagf(z);
-  const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f));
+  const float __pi(static_cast<float>(__DEVICELIB_PI));
   if (__spirv_IsInf(z_real)) {
     if (__spirv_IsNan(z_imag))
       return CMPLXF(z_imag, z_real);
@@ -206,19 +206,25 @@ static inline float __complex__ __devicelib_cacosf(float __complex__ z) {
     return CMPLXF(__pi / 2.0f, -z_imag);
   if (z_real == 0 && (z_imag == 0 || __spirv_IsNan(z_imag)))
     return CMPLXF(__pi / 2.0f, -z_imag);
+  // acos(z) = pi - acos(-z) for real(z) < 0; evaluating on the right
+  // half-plane keeps z + sqrt(z^2-1) from cancelling.  Im(acos) keeps the
+  // sign opposite to imag(z) either way.
+  const bool refl = __spirv_SignBitSet(z_real);
+  float __complex__ zr = refl ? -z : z;
   float __complex__ w =
-      __devicelib_clogf(z + __devicelib_csqrtf(__sqrf(z) - 1.0f));
+      __devicelib_clogf(zr + __devicelib_csqrtf(__sqrf(zr) - 1.0f));
+  float re = __spirv_ocl_fabs(__devicelib_cimagf(w));
+  if (refl)
+    re = __pi - re;
   if (__spirv_SignBitSet(z_imag))
-    return CMPLXF(__spirv_ocl_fabs(__devicelib_cimagf(w)),
-                  __spirv_ocl_fabs(__devicelib_crealf(w)));
-  return CMPLXF(__spirv_ocl_fabs(__devicelib_cimagf(w)),
-                -__spirv_ocl_fabs(__devicelib_crealf(w)));
+    return CMPLXF(re, __spirv_ocl_fabs(__devicelib_crealf(w)));
+  return CMPLXF(re, -__spirv_ocl_fabs(__devicelib_crealf(w)));
 }
 
 static inline float __complex__ __devicelib_casinhf(float __complex__ z) {
   float z_real = __devicelib_crealf(z);
   float z_imag = __devicelib_cimagf(z);
-  const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f));
+  const float __pi(static_cast<float>(__DEVICELIB_PI));
   if (__spirv_IsInf(z_real)) {
     if (__spirv_IsNan(z_imag))
       return z;
@@ -236,8 +242,11 @@ static inline float __complex__ __devicelib_casinhf(float __complex__ z) {
   if (__spirv_IsInf(z_imag))
     return CMPLXF(__spirv_ocl_copysign(z_imag, z_real),
                   __spirv_ocl_copysign(__pi / 2.0f, z_imag));
+  // asinh is odd: evaluate on the right half-plane, where sqrt(z^2+1) ~ +z and
+  // the sum cannot cancel, then take the signs from the original argument.
+  float __complex__ zr = __spirv_SignBitSet(z_real) ? -z : z;
   float __complex__ w =
-      __devicelib_clogf(z + __devicelib_csqrtf(__sqrf(z) + 1.0f));
+      __devicelib_clogf(zr + __devicelib_csqrtf(__sqrf(zr) + 1.0f));
   return CMPLXF(__spirv_ocl_copysign(__devicelib_crealf(w), z_real),
                 __spirv_ocl_copysign(__devicelib_cimagf(w), z_imag));
 }
@@ -251,7 +260,7 @@ static inline float __complex__ __devicelib_casinf(float __complex__ z) {
 static inline float __complex__ __devicelib_cacoshf(float __complex__ z) {
   float z_real = __devicelib_crealf(z);
   float z_imag = __devicelib_cimagf(z);
-  const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f));
+  const float __pi(static_cast<float>(__DEVICELIB_PI));
   if (__spirv_IsInf(z_real)) {
     if (__spirv_IsNan(z_imag))
       return CMPLXF(__spirv_ocl_fabs(z_real), z_imag);
@@ -273,16 +282,23 @@ static inline float __complex__ __devicelib_cacoshf(float __complex__ z) {
   if (__spirv_IsInf(z_imag))
     return CMPLXF(__spirv_ocl_fabs(z_imag),
                   __spirv_ocl_copysign(__pi / 2.0f, z_imag));
+  // acosh(z) = acosh(-z) + i*copysign(pi, imag(z)) for real(z) < 0; evaluating
+  // on the right half-plane keeps z + sqrt(z^2-1) from cancelling.
+  const bool refl = __spirv_SignBitSet(z_real);
+  float __complex__ zr = refl ? -z : z;
   float __complex__ w =
-      __devicelib_clogf(z + __devicelib_csqrtf(__sqrf(z) - 1.0f));
+      __devicelib_clogf(zr + __devicelib_csqrtf(__sqrf(zr) - 1.0f));
+  float im = __spirv_ocl_fabs(__devicelib_cimagf(w));
+  if (refl)
+    im = __pi - im;
   return CMPLXF(__spirv_ocl_copysign(__devicelib_crealf(w), 0.0f),
-                __spirv_ocl_copysign(__devicelib_cimagf(w), z_imag));
+                __spirv_ocl_copysign(im, z_imag));
 }
 
 static inline float __complex__ __devicelib_catanhf(float __complex__ z) {
   float z_real = __devicelib_crealf(z);
   float z_imag = __devicelib_cimagf(z);
-  const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f));
+  const float __pi(static_cast<float>(__DEVICELIB_PI));
   if (__spirv_IsInf(z_imag))
     return CMPLXF(__spirv_ocl_copysign(0.0f, z_real),
                   __spirv_ocl_copysign(__pi / 2.0f, z_imag));

sycl/include/sycl/detail/numbers.hpp

@@ -0,0 +1,32 @@
+//==---------------- numbers.hpp - math constants --------------------------==//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#pragma once
+
+#if __cplusplus >= 202002L && __has_include(<version>)
+#include <version> // defines __cpp_lib_math_constants
+#endif
+#if __cpp_lib_math_constants
+#include <numbers>
+#endif
+
+namespace sycl {
+inline namespace _V1 {
+namespace detail::numbers {
+
+// std::numbers::pi once C++20 is the floor. Held as double rather than
+// pi_v<T>: pi_v<sycl::half> is ill-formed, so callers convert at the use site.
+#if __cpp_lib_math_constants
+inline constexpr double pi = std::numbers::pi;
+#else
+inline constexpr double pi = 3.14159265358979323846;
+#endif
+
+} // namespace detail::numbers
+} // namespace _V1
+} // namespace sycl

sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp

@@ -11,6 +11,7 @@
 #include "common.hpp"
 
 #include <sycl/builtins.hpp>
+#include <sycl/detail/numbers.hpp>
 
 #include <math.h>
 
@@ -357,7 +358,7 @@ template <class _Tp>
 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
     typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
     asinh(const complex<_Tp> &__x) {
-  const _Tp __pi(sycl::atan2(_Tp(+0.), _Tp(-0.)));
+  const _Tp __pi(static_cast<_Tp>(sycl::detail::numbers::pi));
   if (sycl::isinf(__x.real())) {
     if (sycl::isnan(__x.imag()))
       return __x;
@@ -376,7 +377,10 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
   if (sycl::isinf(__x.imag()))
     return complex<_Tp>(sycl::copysign(__x.imag(), __x.real()),
                         sycl::copysign(__pi / _Tp(2), __x.imag()));
-  complex<_Tp> __z = log(__x + sqrt(cplx::detail::__sqr(__x) + _Tp(1)));
+  // asinh is odd: evaluate on the right half-plane, where sqrt(x^2+1) ~ +x and
+  // the sum cannot cancel, then take the signs from the original argument.
+  const complex<_Tp> __xr = sycl::signbit(__x.real()) ? -__x : __x;
+  complex<_Tp> __z = log(__xr + sqrt(cplx::detail::__sqr(__xr) + _Tp(1)));
   return complex<_Tp>(sycl::copysign(__z.real(), __x.real()),
                       sycl::copysign(__z.imag(), __x.imag()));
 }
@@ -387,7 +391,7 @@ template <class _Tp>
 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
     typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
     acosh(const complex<_Tp> &__x) {
-  const _Tp __pi(sycl::atan2(_Tp(+0.), _Tp(-0.)));
+  const _Tp __pi(static_cast<_Tp>(sycl::detail::numbers::pi));
   if (sycl::isinf(__x.real())) {
     if (sycl::isnan(__x.imag()))
       return complex<_Tp>(sycl::fabs(__x.real()), __x.imag());
@@ -411,9 +415,15 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
   if (sycl::isinf(__x.imag()))
     return complex<_Tp>(sycl::fabs(__x.imag()),
                         sycl::copysign(__pi / _Tp(2), __x.imag()));
-  complex<_Tp> __z = log(__x + sqrt(cplx::detail::__sqr(__x) - _Tp(1)));
+  // acosh(x) = acosh(-x) + i*copysign(pi, imag(x)) for real(x) < 0; evaluating
+  // on the right half-plane keeps x + sqrt(x^2-1) from cancelling.
+  const bool __refl = sycl::signbit(__x.real());
+  const complex<_Tp> __xr = __refl ? -__x : __x;
+  complex<_Tp> __z = log(__xr + sqrt(cplx::detail::__sqr(__xr) - _Tp(1)));
+  const _Tp __im =
+      __refl ? __pi - sycl::fabs(__z.imag()) : sycl::fabs(__z.imag());
   return complex<_Tp>(sycl::copysign(__z.real(), _Tp(0)),
-                      sycl::copysign(__z.imag(), __x.imag()));
+                      sycl::copysign(__im, __x.imag()));
 }
 
 // atanh
@@ -422,7 +432,7 @@ template <class _Tp>
 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
     typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
     atanh(const complex<_Tp> &__x) {
-  const _Tp __pi(sycl::atan2(_Tp(+0.), _Tp(-0.)));
+  const _Tp __pi(static_cast<_Tp>(sycl::detail::numbers::pi));
   if (sycl::isinf(__x.imag())) {
     return complex<_Tp>(sycl::copysign(_Tp(0), __x.real()),
                         sycl::copysign(__pi / _Tp(2), __x.imag()));
@@ -526,7 +536,7 @@ template <class _Tp>
 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
     typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
     acos(const complex<_Tp> &__x) {
-  const _Tp __pi(sycl::atan2(_Tp(+0.), _Tp(-0.)));
+  const _Tp __pi(static_cast<_Tp>(sycl::detail::numbers::pi));
   if (sycl::isinf(__x.real())) {
     if (sycl::isnan(__x.imag()))
       return complex<_Tp>(__x.imag(), __x.real());
@@ -550,10 +560,17 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
     return complex<_Tp>(__pi / _Tp(2), -__x.imag());
   if (__x.real() == 0 && (__x.imag() == 0 || sycl::isnan(__x.imag())))
     return complex<_Tp>(__pi / _Tp(2), -__x.imag());
-  complex<_Tp> __z = log(__x + sqrt(cplx::detail::__sqr(__x) - _Tp(1)));
+  // acos(x) = pi - acos(-x) for real(x) < 0; evaluating on the right
+  // half-plane keeps x + sqrt(x^2-1) from cancelling.  Im(acos) keeps the
+  // sign opposite to imag(x) either way.
+  const bool __refl = sycl::signbit(__x.real());
+  const complex<_Tp> __xr = __refl ? -__x : __x;
+  complex<_Tp> __z = log(__xr + sqrt(cplx::detail::__sqr(__xr) - _Tp(1)));
+  const _Tp __re =
+      __refl ? __pi - sycl::fabs(__z.imag()) : sycl::fabs(__z.imag());
   if (sycl::signbit(__x.imag()))
-    return complex<_Tp>(sycl::fabs(__z.imag()), sycl::fabs(__z.real()));
-  return complex<_Tp>(sycl::fabs(__z.imag()), -sycl::fabs(__z.real()));
+    return complex<_Tp>(__re, sycl::fabs(__z.real()));
+  return complex<_Tp>(__re, -sycl::fabs(__z.real()));
 }
 
 // atan