From dcb335db9383a9c3b5e7c09851972697c1f9a199 Mon Sep 17 00:00:00 2001 From: Chris Perkins Date: Mon, 5 Oct 2026 11:59:47 -0700 Subject: [PATCH 1/4] improve complex trig and pi ulp Signed-off-by: Chris Perkins --- libdevice/device_complex.h | 12 +++ libdevice/fallback-complex-fp64.hpp | 40 +++++++--- libdevice/fallback-complex.hpp | 40 +++++++--- .../complex/detail/complex_math.hpp | 51 ++++++++++--- .../Complex/sycl_complex_math_test_cases.hpp | 34 ++++++++- .../std_complex_inverse_trig_large_args.cpp | 76 +++++++++++++++++++ 6 files changed, 217 insertions(+), 36 deletions(-) create mode 100644 sycl/test-e2e/DeviceLib/std_complex_inverse_trig_large_args.cpp diff --git a/libdevice/device_complex.h b/libdevice/device_complex.h index a7be1421ee575..63ca6079453e2 100644 --- a/libdevice/device_complex.h +++ b/libdevice/device_complex.h @@ -10,6 +10,18 @@ #include "device.h" +#if __cplusplus >= 202002L && __has_include() +#include // 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 +#define __DEVICELIB_PI std::numbers::pi_v +#else +#define __DEVICELIB_PI 3.14159265358979323846 +#endif + #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..3083a7db97e31 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)); + 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)); diff --git a/libdevice/fallback-complex.hpp b/libdevice/fallback-complex.hpp index 6b9961aa94c05..9a186981c763c 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)); + const float __pi(static_cast(__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(__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(__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(__DEVICELIB_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/ext/oneapi/experimental/complex/detail/complex_math.hpp b/sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp index c9557f83d2093..42eadfac947a3 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 @@ -14,6 +14,13 @@ #include +#if __cplusplus >= 202002L && __has_include() +#include // defines __cpp_lib_math_constants +#endif +#if __cpp_lib_math_constants +#include +#endif + namespace sycl { inline namespace _V1 { @@ -340,6 +347,14 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY } namespace cplx::detail { +// Correctly rounded pi; the atan2(+0, -0) idiom from libc++ is 1 ulp low for +// float on Intel GPUs. Kept as double: pi_v is ill-formed. +#if __cpp_lib_math_constants +inline constexpr double __pi_d = std::numbers::pi_v; +#else +inline constexpr double __pi_d = 3.14159265358979323846; +#endif + // __sqr, computes pow(x, 2) template @@ -357,7 +372,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>(cplx::detail::__pi_d)); if (sycl::isinf(__x.real())) { if (sycl::isnan(__x.imag())) return __x; @@ -376,7 +391,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 +405,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>(cplx::detail::__pi_d)); if (sycl::isinf(__x.real())) { if (sycl::isnan(__x.imag())) return complex<_Tp>(sycl::fabs(__x.real()), __x.imag()); @@ -411,9 +429,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 +446,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>(cplx::detail::__pi_d)); if (sycl::isinf(__x.imag())) { return complex<_Tp>(sycl::copysign(_Tp(0), __x.real()), sycl::copysign(__pi / _Tp(2), __x.imag())); @@ -526,7 +550,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>(cplx::detail::__pi_d)); if (sycl::isinf(__x.real())) { if (sycl::isnan(__x.imag())) return complex<_Tp>(__x.imag(), __x.real()); @@ -550,10 +574,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 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