diff --git a/libdevice/cmake/modules/SYCLLibdevice.cmake b/libdevice/cmake/modules/SYCLLibdevice.cmake index 5bb9689a74bf2..f7d996923dd85 100644 --- a/libdevice/cmake/modules/SYCLLibdevice.cmake +++ b/libdevice/cmake/modules/SYCLLibdevice.cmake @@ -327,7 +327,8 @@ set(sycl-compiler_deps set(crt_obj_deps wrapper.h device.h spirv_vars.h fallback-cstring.hpp ${sycl-compiler_deps}) set(cmath_obj_deps device_math.h device.h fallback-cmath.hpp fallback-cmath-fp64.hpp - device_complex.h fallback-complex.hpp fallback-complex-fp64.hpp complex_wrapper.hpp ${sycl-compiler_deps}) + device_complex.h fallback-complex.hpp fallback-complex-fp64.hpp complex_wrapper.hpp + sycl-headers ${sycl-compiler_deps}) set(imf_obj_deps device_imf.hpp imf_half.hpp imf_bf16.hpp imf_rounding_op.hpp imf_impl_utils.hpp device.h ${sycl-compiler_deps}) set(itt_obj_deps device_itt.h spirv_vars.h device.h ${sycl-compiler_deps}) set(bfloat16_obj_deps sycl-headers ${sycl-compiler_deps}) diff --git a/libdevice/device_complex.h b/libdevice/device_complex.h index a7be1421ee575..70a5a4a2bdff2 100644 --- a/libdevice/device_complex.h +++ b/libdevice/device_complex.h @@ -10,6 +10,8 @@ #include "device.h" +#include + #if defined(__SPIR__) || defined(__SPIRV__) // TODO: This needs to be more robust. diff --git a/libdevice/fallback-complex-fp64.hpp b/libdevice/fallback-complex-fp64.hpp index e2d80b30dd3a1..910f621d5f50c 100644 --- a/libdevice/fallback-complex-fp64.hpp +++ b/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)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return CMPLX(z_imag, z_real); @@ -213,19 +213,23 @@ 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); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use acos(z) = pi - acos(-z). + 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)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return z; @@ -243,8 +247,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)); + // For real(z) < 0, z + sqrt(z^2+1) cancels; use asinh(z) = -asinh(-z), + // with the copysign calls below restoring the signs. + 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 +265,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)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return CMPLX(__spirv_ocl_fabs(z_real), z_imag); @@ -280,16 +287,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)); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use + // acosh(z) = acosh(-z) + i*copysign(pi, imag(z)). + 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)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_imag)) return CMPLX(__spirv_ocl_copysign(0.0, z_real), __spirv_ocl_copysign(__pi / 2.0, z_imag)); diff --git a/libdevice/fallback-complex.hpp b/libdevice/fallback-complex.hpp index 6b9961aa94c05..74a7280e2e42b 100644 --- a/libdevice/fallback-complex.hpp +++ b/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)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return CMPLXF(z_imag, z_real); @@ -206,19 +206,23 @@ 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); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use acos(z) = pi - acos(-z). + 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)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return z; @@ -236,8 +240,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)); + // For real(z) < 0, z + sqrt(z^2+1) cancels; use asinh(z) = -asinh(-z), + // with the copysign calls below restoring the signs. + 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 +258,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)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return CMPLXF(__spirv_ocl_fabs(z_real), z_imag); @@ -273,16 +280,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)); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use + // acosh(z) = acosh(-z) + i*copysign(pi, imag(z)). + 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)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_imag)) return CMPLXF(__spirv_ocl_copysign(0.0f, z_real), __spirv_ocl_copysign(__pi / 2.0f, z_imag)); diff --git a/sycl/include/sycl/detail/numbers.hpp b/sycl/include/sycl/detail/numbers.hpp new file mode 100644 index 0000000000000..269629c7b44c3 --- /dev/null +++ b/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() +#include // defines __cpp_lib_math_constants +#endif +#if __cpp_lib_math_constants +#include +#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: pi_v 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 diff --git a/sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp b/sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp index c9557f83d2093..857474bc9e71a 100644 --- a/sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp +++ b/sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp @@ -11,6 +11,7 @@ #include "common.hpp" #include +#include #include @@ -357,7 +358,7 @@ template __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY typename std::enable_if_t::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))); + // For real(x) < 0, x + sqrt(x^2+1) cancels; use asinh(x) = -asinh(-x), + // with the copysign calls below restoring the signs. + 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 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY typename std::enable_if_t::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))); + // For real(x) < 0, x + sqrt(x^2-1) cancels; use + // acosh(x) = acosh(-x) + i*copysign(pi, imag(x)). + 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 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY typename std::enable_if_t::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 __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY typename std::enable_if_t::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,15 @@ __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))); + // For real(x) < 0, x + sqrt(x^2-1) cancels; use acos(x) = pi - acos(-x). + 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 diff --git a/sycl/include/sycl/stl_wrappers/__sycl_complex_impl.hpp b/sycl/include/sycl/stl_wrappers/__sycl_complex_impl.hpp index 83d23a99b7477..00b910c0e50b4 100644 --- a/sycl/include/sycl/stl_wrappers/__sycl_complex_impl.hpp +++ b/sycl/include/sycl/stl_wrappers/__sycl_complex_impl.hpp @@ -14,6 +14,7 @@ #ifdef __SYCL_DEVICE_ONLY__ #include +#include // The 'sycl_device_only' attribute enables device-side overloading. #define __SYCL_DEVICE __attribute__((sycl_device_only, always_inline)) @@ -603,7 +604,7 @@ __SYCL_DEVICE_C float __complex__ cacosf(float __complex__ z) { float z_real = crealf(z); float z_imag = cimagf(z); - const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return __SYCL_CMPLXF(z_imag, z_real); @@ -625,18 +626,22 @@ float __complex__ cacosf(float __complex__ z) { return __SYCL_CMPLXF(__pi / 2.0f, -z_imag); if (z_real == 0 && (z_imag == 0 || __spirv_IsNan(z_imag))) return __SYCL_CMPLXF(__pi / 2.0f, -z_imag); - float __complex__ w = clogf(z + csqrtf(__sqrf(z) - 1.0f)); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use acos(z) = pi - acos(-z). + const bool refl = __spirv_SignBitSet(z_real); + float __complex__ zr = refl ? -z : z; + float __complex__ w = clogf(zr + csqrtf(__sqrf(zr) - 1.0f)); + float re = __spirv_ocl_fabs(cimagf(w)); + if (refl) + re = __pi - re; if (__spirv_SignBitSet(z_imag)) - return __SYCL_CMPLXF(__spirv_ocl_fabs(cimagf(w)), - __spirv_ocl_fabs(crealf(w))); - return __SYCL_CMPLXF(__spirv_ocl_fabs(cimagf(w)), - -__spirv_ocl_fabs(crealf(w))); + return __SYCL_CMPLXF(re, __spirv_ocl_fabs(crealf(w))); + return __SYCL_CMPLXF(re, -__spirv_ocl_fabs(crealf(w))); } __SYCL_DEVICE_C double __complex__ cacos(double __complex__ z) { double z_real = creal(z); double z_imag = cimag(z); - const double __pi(__spirv_ocl_atan2(+0.0, -0.0)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return __SYCL_CMPLX(z_imag, z_real); @@ -658,17 +663,23 @@ double __complex__ cacos(double __complex__ z) { return __SYCL_CMPLX(__pi / 2.0, -z_imag); if (z_real == 0 && (z_imag == 0 || __spirv_IsNan(z_imag))) return __SYCL_CMPLX(__pi / 2.0, -z_imag); - double __complex__ w = clog(z + csqrt(__sqr(z) - 1.0)); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use acos(z) = pi - acos(-z). + const bool refl = __spirv_SignBitSet(z_real); + double __complex__ zr = refl ? -z : z; + double __complex__ w = clog(zr + csqrt(__sqr(zr) - 1.0)); + double re = __spirv_ocl_fabs(cimag(w)); + if (refl) + re = __pi - re; if (__spirv_SignBitSet(z_imag)) - return __SYCL_CMPLX(__spirv_ocl_fabs(cimag(w)), __spirv_ocl_fabs(creal(w))); - return __SYCL_CMPLX(__spirv_ocl_fabs(cimag(w)), -__spirv_ocl_fabs(creal(w))); + return __SYCL_CMPLX(re, __spirv_ocl_fabs(creal(w))); + return __SYCL_CMPLX(re, -__spirv_ocl_fabs(creal(w))); } __SYCL_DEVICE_C float __complex__ casinhf(float __complex__ z) { float z_real = crealf(z); float z_imag = cimagf(z); - const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return z; @@ -686,7 +697,10 @@ float __complex__ casinhf(float __complex__ z) { if (__spirv_IsInf(z_imag)) return __SYCL_CMPLXF(__spirv_ocl_copysign(z_imag, z_real), __spirv_ocl_copysign(__pi / 2.0f, z_imag)); - float __complex__ w = clogf(z + csqrtf(__sqrf(z) + 1.0f)); + // For real(z) < 0, z + sqrt(z^2+1) cancels; use asinh(z) = -asinh(-z), + // with the copysign calls below restoring the signs. + float __complex__ zr = __spirv_SignBitSet(z_real) ? -z : z; + float __complex__ w = clogf(zr + csqrtf(__sqrf(zr) + 1.0f)); return __SYCL_CMPLXF(__spirv_ocl_copysign(crealf(w), z_real), __spirv_ocl_copysign(cimagf(w), z_imag)); } @@ -694,7 +708,7 @@ __SYCL_DEVICE_C double __complex__ casinh(double __complex__ z) { double z_real = creal(z); double z_imag = cimag(z); - const double __pi(__spirv_ocl_atan2(+0.0, -0.0)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return z; @@ -712,7 +726,10 @@ double __complex__ casinh(double __complex__ z) { if (__spirv_IsInf(z_imag)) return __SYCL_CMPLX(__spirv_ocl_copysign(z_imag, z_real), __spirv_ocl_copysign(__pi / 2.0, z_imag)); - double __complex__ w = clog(z + csqrt(__sqr(z) + 1.0)); + // For real(z) < 0, z + sqrt(z^2+1) cancels; use asinh(z) = -asinh(-z), + // with the copysign calls below restoring the signs. + double __complex__ zr = __spirv_SignBitSet(z_real) ? -z : z; + double __complex__ w = clog(zr + csqrt(__sqr(zr) + 1.0)); return __SYCL_CMPLX(__spirv_ocl_copysign(creal(w), z_real), __spirv_ocl_copysign(cimag(w), z_imag)); } @@ -732,7 +749,7 @@ __SYCL_DEVICE_C float __complex__ cacoshf(float __complex__ z) { float z_real = crealf(z); float z_imag = cimagf(z); - const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return __SYCL_CMPLXF(__spirv_ocl_fabs(z_real), z_imag); @@ -756,15 +773,22 @@ float __complex__ cacoshf(float __complex__ z) { if (__spirv_IsInf(z_imag)) return __SYCL_CMPLXF(__spirv_ocl_fabs(z_imag), __spirv_ocl_copysign(__pi / 2.0f, z_imag)); - float __complex__ w = clogf(z + csqrtf(__sqrf(z) - 1.0f)); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use + // acosh(z) = acosh(-z) + i*copysign(pi, imag(z)). + const bool refl = __spirv_SignBitSet(z_real); + float __complex__ zr = refl ? -z : z; + float __complex__ w = clogf(zr + csqrtf(__sqrf(zr) - 1.0f)); + float im = __spirv_ocl_fabs(cimagf(w)); + if (refl) + im = __pi - im; return __SYCL_CMPLXF(__spirv_ocl_copysign(crealf(w), 0.0f), - __spirv_ocl_copysign(cimagf(w), z_imag)); + __spirv_ocl_copysign(im, z_imag)); } __SYCL_DEVICE_C double __complex__ cacosh(double __complex__ z) { double z_real = creal(z); double z_imag = cimag(z); - const double __pi(__spirv_ocl_atan2(+0.0, -0.0)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_real)) { if (__spirv_IsNan(z_imag)) return __SYCL_CMPLX(__spirv_ocl_fabs(z_real), z_imag); @@ -787,16 +811,23 @@ double __complex__ cacosh(double __complex__ z) { if (__spirv_IsInf(z_imag)) return __SYCL_CMPLX(__spirv_ocl_fabs(z_imag), __spirv_ocl_copysign(__pi / 2.0, z_imag)); - double __complex__ w = clog(z + csqrt(__sqr(z) - 1.0)); + // For real(z) < 0, z + sqrt(z^2-1) cancels; use + // acosh(z) = acosh(-z) + i*copysign(pi, imag(z)). + const bool refl = __spirv_SignBitSet(z_real); + double __complex__ zr = refl ? -z : z; + double __complex__ w = clog(zr + csqrt(__sqr(zr) - 1.0)); + double im = __spirv_ocl_fabs(cimag(w)); + if (refl) + im = __pi - im; return __SYCL_CMPLX(__spirv_ocl_copysign(creal(w), 0.0), - __spirv_ocl_copysign(cimag(w), z_imag)); + __spirv_ocl_copysign(im, z_imag)); } __SYCL_DEVICE_C float __complex__ catanhf(float __complex__ z) { float z_real = crealf(z); float z_imag = cimagf(z); - const float __pi(__spirv_ocl_atan2(+0.0f, -0.0f)); + constexpr float __pi = static_cast(sycl::detail::numbers::pi); if (__spirv_IsInf(z_imag)) return __SYCL_CMPLXF(__spirv_ocl_copysign(0.0f, z_real), __spirv_ocl_copysign(__pi / 2.0f, z_imag)); @@ -826,7 +857,7 @@ __SYCL_DEVICE_C double __complex__ catanh(double __complex__ z) { double z_real = creal(z); double z_imag = cimag(z); - const double __pi(__spirv_ocl_atan2(+0.0, -0.0)); + constexpr double __pi = sycl::detail::numbers::pi; if (__spirv_IsInf(z_imag)) return __SYCL_CMPLX(__spirv_ocl_copysign(0.0, z_real), __spirv_ocl_copysign(__pi / 2.0, z_imag)); diff --git a/sycl/test-e2e/Complex/sycl_complex_math_test_cases.hpp b/sycl/test-e2e/Complex/sycl_complex_math_test_cases.hpp index d2616c5e9be92..49d9bc45ecc52 100644 --- a/sycl/test-e2e/Complex/sycl_complex_math_test_cases.hpp +++ b/sycl/test-e2e/Complex/sycl_complex_math_test_cases.hpp @@ -103,8 +103,7 @@ struct deci_test_cases { } }; -template