Skip to content

[SYCL] Use numeric_limits instead of NAN And INFINITY macros #19852

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Sep 1, 2025
Merged

[SYCL] Use numeric_limits instead of NAN And INFINITY macros #19852

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fccc345
[SYCL] Use numeric_limits instead of NAN And INFINITY macros
steffenlarsen Aug 22, 2025
9ce6cb7
Fix missed NAN
steffenlarsen Aug 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion sycl/include/sycl/ext/oneapi/experimental/bfloat16_math.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ constexpr int num_elements_v = sycl::detail::num_elements<T>::value;

/******************* isnan ********************/

// According to bfloat16 format, NAN value's exponent field is 0xFF and
// According to bfloat16 format, NaN value's exponent field is 0xFF and
// significand has non-zero bits.
template <typename T>
std::enable_if_t<std::is_same_v<T, bfloat16>, bool> isnan(T x) {
Expand Down
33 changes: 21 additions & 12 deletions sycl/include/sycl/ext/oneapi/experimental/complex/detail/complex_math.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,8 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
proj(const complex<_Tp> &__c) {
complex<_Tp> __r = __c;
if (sycl::isinf(__c.real()) || sycl::isinf(__c.imag()))
__r = complex<_Tp>(INFINITY, sycl::copysign(_Tp(0), __c.imag()));
__r = complex<_Tp>(std::numeric_limits<float>::infinity(),
sycl::copysign(_Tp(0), __c.imag()));
return __r;
}

Expand All @@ -214,16 +215,18 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
polar(const _Tp &__rho, const _Tp &__theta = _Tp()) {
if (sycl::isnan(__rho) || sycl::signbit(__rho))
return complex<_Tp>(_Tp(NAN), _Tp(NAN));
return complex<_Tp>(_Tp(std::numeric_limits<float>::quiet_NaN()),
_Tp(std::numeric_limits<float>::quiet_NaN()));
if (sycl::isnan(__theta)) {
if (sycl::isinf(__rho))
return complex<_Tp>(__rho, __theta);
return complex<_Tp>(__theta, __theta);
}
if (sycl::isinf(__theta)) {
if (sycl::isinf(__rho))
return complex<_Tp>(__rho, _Tp(NAN));
return complex<_Tp>(_Tp(NAN), _Tp(NAN));
return complex<_Tp>(__rho, _Tp(std::numeric_limits<float>::quiet_NaN()));
return complex<_Tp>(_Tp(std::numeric_limits<float>::quiet_NaN()),
_Tp(std::numeric_limits<float>::quiet_NaN()));
}
_Tp __x = __rho * sycl::cos(__theta);
if (sycl::isnan(__x))
Expand Down Expand Up @@ -259,7 +262,8 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
sqrt(const complex<_Tp> &__x) {
if (sycl::isinf(__x.imag()))
return complex<_Tp>(_Tp(INFINITY), __x.imag());
return complex<_Tp>(_Tp(std::numeric_limits<float>::infinity()),
__x.imag());
if (sycl::isinf(__x.real())) {
if (__x.real() > _Tp(0))
return complex<_Tp>(__x.real(), sycl::isnan(__x.imag())
Expand Down Expand Up @@ -288,7 +292,7 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
__i = _Tp(1);
} else if (__i == 0 || !sycl::isfinite(__i)) {
if (sycl::isinf(__i))
__i = _Tp(NAN);
__i = _Tp(std::numeric_limits<float>::quiet_NaN());
return complex<_Tp>(__x.real(), __i);
}
}
Expand Down Expand Up @@ -436,8 +440,9 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
sycl::copysign(__pi / _Tp(2), __x.imag()));
}
if (sycl::fabs(__x.real()) == _Tp(1) && __x.imag() == _Tp(0)) {
return complex<_Tp>(sycl::copysign(_Tp(INFINITY), __x.real()),
sycl::copysign(_Tp(0), __x.imag()));
return complex<_Tp>(
sycl::copysign(_Tp(std::numeric_limits<float>::infinity()), __x.real()),
sycl::copysign(_Tp(0), __x.imag()));
}
complex<_Tp> __z = log((_Tp(1) + __x) / (_Tp(1) - __x)) / _Tp(2);
return complex<_Tp>(sycl::copysign(__z.real(), __x.real()),
Expand All @@ -451,9 +456,11 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
sinh(const complex<_Tp> &__x) {
if (sycl::isinf(__x.real()) && !sycl::isfinite(__x.imag()))
return complex<_Tp>(__x.real(), _Tp(NAN));
return complex<_Tp>(__x.real(),
_Tp(std::numeric_limits<float>::quiet_NaN()));
if (__x.real() == 0 && !sycl::isfinite(__x.imag()))
return complex<_Tp>(__x.real(), _Tp(NAN));
return complex<_Tp>(__x.real(),
_Tp(std::numeric_limits<float>::quiet_NaN()));
if (__x.imag() == 0 && !sycl::isfinite(__x.real()))
return __x;
return complex<_Tp>(sycl::sinh(__x.real()) * sycl::cos(__x.imag()),
Expand All @@ -467,9 +474,11 @@ __DPCPP_SYCL_EXTERNAL _SYCL_EXT_CPLX_INLINE_VISIBILITY
typename std::enable_if_t<is_genfloat<_Tp>::value, complex<_Tp>>
cosh(const complex<_Tp> &__x) {
if (sycl::isinf(__x.real()) && !sycl::isfinite(__x.imag()))
return complex<_Tp>(sycl::fabs(__x.real()), _Tp(NAN));
return complex<_Tp>(sycl::fabs(__x.real()),
_Tp(std::numeric_limits<float>::quiet_NaN()));
if (__x.real() == 0 && !sycl::isfinite(__x.imag()))
return complex<_Tp>(_Tp(NAN), __x.real());
return complex<_Tp>(_Tp(std::numeric_limits<float>::quiet_NaN()),
__x.real());
if (__x.real() == 0 && __x.imag() == 0)
return complex<_Tp>(_Tp(1), __x.imag());
if (__x.imag() == 0 && !sycl::isfinite(__x.real()))
Expand Down
98 changes: 60 additions & 38 deletions sycl/include/sycl/stl_wrappers/__sycl_complex_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@
// This header defines device-side overloads of <complex> functions.

#ifdef __SYCL_DEVICE_ONLY__
#include <cmath> // For INFINITY.

#include <limits>

// The 'sycl_device_only' attribute enables device-side overloading.
#define __SYCL_DEVICE __attribute__((sycl_device_only, always_inline))
Expand Down Expand Up @@ -82,8 +83,9 @@ float __complex__ __mulsc3(float __a, float __b, float __c, float __d) {
__recalc = 1.0f;
}
if (__recalc) {
z = __SYCL_CMPLXF((INFINITY * (__a * __c - __b * __d)),
(INFINITY * (__a * __d + __b * __c)));
z = __SYCL_CMPLXF(
(std::numeric_limits<float>::infinity() * (__a * __c - __b * __d)),
(std::numeric_limits<float>::infinity() * (__a * __d + __b * __c)));
}
}
return z;
Expand Down Expand Up @@ -131,8 +133,9 @@ double __complex__ __muldc3(double __a, double __b, double __c, double __d) {
__recalc = 1.0;
}
if (__recalc) {
z = __SYCL_CMPLX((INFINITY * (__a * __c - __b * __d)),
(INFINITY * (__a * __d + __b * __c)));
z = __SYCL_CMPLX(
(std::numeric_limits<float>::infinity() * (__a * __c - __b * __d)),
Copy link
Contributor

@cperkinsintel cperkinsintel Aug 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't this be std::numeric_limits<double> (instead of float)?

Copy link
Contributor Author

@steffenlarsen steffenlarsen Aug 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The only cases I changed the type were the ones with an immediate cast after them. Though I could change it here, I am worried that there could be some unexpected changes in rounding as a result of that.

(std::numeric_limits<float>::infinity() * (__a * __d + __b * __c)));
}
}
return z;
Expand Down Expand Up @@ -161,15 +164,19 @@ float __complex__ __divsc3(float __a, float __b, float __c, float __d) {
z = __SYCL_CMPLXF(z_real, z_imag);
if (__spirv_IsNan(z_real) && __spirv_IsNan(z_imag)) {
if ((__denom == 0.0f) && (!__spirv_IsNan(__a) || !__spirv_IsNan(__b))) {
z_real = __spirv_ocl_copysign(INFINITY, __c) * __a;
z_imag = __spirv_ocl_copysign(INFINITY, __c) * __b;
z_real =
__spirv_ocl_copysign(std::numeric_limits<float>::infinity(), __c) *
__a;
z_imag =
__spirv_ocl_copysign(std::numeric_limits<float>::infinity(), __c) *
__b;
z = __SYCL_CMPLXF(z_real, z_imag);
} else if ((__spirv_IsInf(__a) || __spirv_IsInf(__b)) &&
__spirv_IsFinite(__c) && __spirv_IsFinite(__d)) {
__a = __spirv_ocl_copysign(__spirv_IsInf(__a) ? 1.0f : 0.0f, __a);
__b = __spirv_ocl_copysign(__spirv_IsInf(__b) ? 1.0f : 0.0f, __b);
z_real = INFINITY * (__a * __c + __b * __d);
z_imag = INFINITY * (__b * __c - __a * __d);
z_real = std::numeric_limits<float>::infinity() * (__a * __c + __b * __d);
z_imag = std::numeric_limits<float>::infinity() * (__b * __c - __a * __d);
z = __SYCL_CMPLXF(z_real, z_imag);
} else if (__spirv_IsInf(__logbw) && __logbw > 0.0f &&
__spirv_IsFinite(__a) && __spirv_IsFinite(__b)) {
Expand Down Expand Up @@ -203,15 +210,19 @@ double __complex__ __divdc3(double __a, double __b, double __c, double __d) {
z = __SYCL_CMPLX(z_real, z_imag);
if (__spirv_IsNan(z_real) && __spirv_IsNan(z_imag)) {
if ((__denom == 0.0) && (!__spirv_IsNan(__a) || !__spirv_IsNan(__b))) {
z_real = __spirv_ocl_copysign((double)INFINITY, __c) * __a;
z_imag = __spirv_ocl_copysign((double)INFINITY, __c) * __b;
z_real =
__spirv_ocl_copysign(std::numeric_limits<double>::infinity(), __c) *
__a;
z_imag =
__spirv_ocl_copysign(std::numeric_limits<double>::infinity(), __c) *
__b;
z = __SYCL_CMPLX(z_real, z_imag);
} else if ((__spirv_IsInf(__a) || __spirv_IsInf(__b)) &&
__spirv_IsFinite(__c) && __spirv_IsFinite(__d)) {
__a = __spirv_ocl_copysign(__spirv_IsInf(__a) ? 1.0 : 0.0, __a);
__b = __spirv_ocl_copysign(__spirv_IsInf(__b) ? 1.0 : 0.0, __b);
z_real = INFINITY * (__a * __c + __b * __d);
z_imag = INFINITY * (__b * __c - __a * __d);
z_real = std::numeric_limits<float>::infinity() * (__a * __c + __b * __d);
z_imag = std::numeric_limits<float>::infinity() * (__b * __c - __a * __d);
Copy link
Contributor

@cperkinsintel cperkinsintel Aug 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ibid. double instead float ?

z = __SYCL_CMPLX(z_real, z_imag);
} else if (__spirv_IsInf(__logbw) && __logbw > 0.0 &&
__spirv_IsFinite(__a) && __spirv_IsFinite(__b)) {
Expand Down Expand Up @@ -247,14 +258,16 @@ __SYCL_DEVICE_C
float __complex__ cprojf(float __complex__ z) {
float __complex__ r = z;
if (__spirv_IsInf(crealf(z)) || __spirv_IsInf(cimagf(z)))
r = __SYCL_CMPLXF(INFINITY, __spirv_ocl_copysign(0.0f, cimagf(z)));
r = __SYCL_CMPLXF(std::numeric_limits<float>::infinity(),
__spirv_ocl_copysign(0.0f, cimagf(z)));
return r;
}
__SYCL_DEVICE_C
double __complex__ cproj(double __complex__ z) {
double __complex__ r = z;
if (__spirv_IsInf(creal(z)) || __spirv_IsInf(cimag(z)))
r = __SYCL_CMPLX(INFINITY, __spirv_ocl_copysign(0.0, cimag(z)));
r = __SYCL_CMPLX(std::numeric_limits<float>::infinity(),
Copy link
Contributor

@cperkinsintel cperkinsintel Aug 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

double again?

__spirv_ocl_copysign(0.0, cimag(z)));
return r;
}

Expand All @@ -275,7 +288,7 @@ float __complex__ cexpf(float __complex__ z) {
return z;
} else if (z_imag == 0.f || !__spirv_IsFinite(z_imag)) {
if (__spirv_IsInf(z_imag))
return __SYCL_CMPLXF(z_real, NAN);
return __SYCL_CMPLXF(z_real, std::numeric_limits<float>::quiet_NaN());
}
}

Expand All @@ -293,16 +306,18 @@ double __complex__ cexp(double __complex__ z) {
z_imag = 1.0;
} else if (z_imag == 0.0 || !__spirv_IsFinite(z_imag)) {
if (__spirv_IsInf(z_imag))
z_imag = NAN;
z_imag = std::numeric_limits<float>::quiet_NaN();
return __SYCL_CMPLX(z_real, z_imag);
}
} else if (__spirv_IsNan(z_real)) {
if (z_imag == 0.0)
return z;
return __SYCL_CMPLX(NAN, NAN);
return __SYCL_CMPLX(std::numeric_limits<float>::quiet_NaN(),
std::numeric_limits<float>::quiet_NaN());
} else if (__spirv_IsFinite(z_real) &&
(__spirv_IsNan(z_imag) || __spirv_IsInf(z_imag))) {
return __SYCL_CMPLX(NAN, NAN);
return __SYCL_CMPLX(std::numeric_limits<float>::quiet_NaN(),
std::numeric_limits<float>::quiet_NaN());
}
double __e = __spirv_ocl_exp(z_real);
double ret_real = __e * __spirv_ocl_cos(z_imag);
Expand Down Expand Up @@ -354,16 +369,18 @@ double __complex__ cpow(double __complex__ x, double __complex__ y) {
__SYCL_DEVICE_C
float __complex__ cpolarf(float rho, float theta) {
if (__spirv_IsNan(rho) || __spirv_SignBitSet(rho))
return __SYCL_CMPLXF(NAN, NAN);
return __SYCL_CMPLXF(std::numeric_limits<float>::quiet_NaN(),
std::numeric_limits<float>::quiet_NaN());
if (__spirv_IsNan(theta)) {
if (__spirv_IsInf(rho))
return __SYCL_CMPLXF(rho, theta);
return __SYCL_CMPLXF(theta, theta);
}
if (__spirv_IsInf(theta)) {
if (__spirv_IsInf(rho))
return __SYCL_CMPLXF(rho, NAN);
return __SYCL_CMPLXF(NAN, NAN);
return __SYCL_CMPLXF(rho, std::numeric_limits<float>::quiet_NaN());
return __SYCL_CMPLXF(std::numeric_limits<float>::quiet_NaN(),
std::numeric_limits<float>::quiet_NaN());
}
float x = rho * __spirv_ocl_cos(theta);
if (__spirv_IsNan(x))
Expand All @@ -376,16 +393,18 @@ float __complex__ cpolarf(float rho, float theta) {
__SYCL_DEVICE_C
double __complex__ cpolar(double rho, double theta) {
if (__spirv_IsNan(rho) || __spirv_SignBitSet(rho))
return __SYCL_CMPLX(NAN, NAN);
return __SYCL_CMPLX(std::numeric_limits<float>::quiet_NaN(),
std::numeric_limits<float>::quiet_NaN());
if (__spirv_IsNan(theta)) {
if (__spirv_IsInf(rho))
return __SYCL_CMPLX(rho, theta);
return __SYCL_CMPLX(theta, theta);
}
if (__spirv_IsInf(theta)) {
if (__spirv_IsInf(rho))
return __SYCL_CMPLX(rho, NAN);
return __SYCL_CMPLX(NAN, NAN);
return __SYCL_CMPLX(rho, std::numeric_limits<float>::quiet_NaN());
return __SYCL_CMPLX(std::numeric_limits<float>::quiet_NaN(),
std::numeric_limits<float>::quiet_NaN());
}
double x = rho * __spirv_ocl_cos(theta);
if (__spirv_IsNan(x))
Expand All @@ -401,7 +420,7 @@ float __complex__ csqrtf(float __complex__ z) {
float z_real = crealf(z);
float z_imag = cimagf(z);
if (__spirv_IsInf(z_imag))
return __SYCL_CMPLXF(INFINITY, z_imag);
return __SYCL_CMPLXF(std::numeric_limits<float>::infinity(), z_imag);
if (__spirv_IsInf(z_real)) {
if (z_real > 0.0f)
return __SYCL_CMPLXF(z_real, __spirv_IsNan(z_imag)
Expand All @@ -417,7 +436,7 @@ double __complex__ csqrt(double __complex__ z) {
double z_real = creal(z);
double z_imag = cimag(z);
if (__spirv_IsInf(z_imag))
return __SYCL_CMPLX(INFINITY, z_imag);
return __SYCL_CMPLX(std::numeric_limits<float>::infinity(), z_imag);
if (__spirv_IsInf(z_real)) {
if (z_real > 0.0)
return __SYCL_CMPLX(z_real, __spirv_IsNan(z_imag)
Expand All @@ -434,9 +453,9 @@ float __complex__ csinhf(float __complex__ z) {
float z_real = crealf(z);
float z_imag = cimagf(z);
if (__spirv_IsInf(z_real) && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLXF(z_real, NAN);
return __SYCL_CMPLXF(z_real, std::numeric_limits<float>::quiet_NaN());
if (z_real == 0 && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLXF(z_real, NAN);
return __SYCL_CMPLXF(z_real, std::numeric_limits<float>::quiet_NaN());
if (z_imag == 0 && !__spirv_IsFinite(z_real))
return z;
return __SYCL_CMPLXF(__spirv_ocl_sinh(z_real) * __spirv_ocl_cos(z_imag),
Expand All @@ -447,9 +466,9 @@ double __complex__ csinh(double __complex__ z) {
double z_real = creal(z);
double z_imag = cimag(z);
if (__spirv_IsInf(z_real) && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLX(z_real, NAN);
return __SYCL_CMPLX(z_real, std::numeric_limits<float>::quiet_NaN());
if (z_real == 0 && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLX(z_real, NAN);
return __SYCL_CMPLX(z_real, std::numeric_limits<float>::quiet_NaN());
if (z_imag == 0 && !__spirv_IsFinite(z_real))
return z;
return __SYCL_CMPLX(__spirv_ocl_sinh(z_real) * __spirv_ocl_cos(z_imag),
Expand All @@ -461,9 +480,10 @@ float __complex__ ccoshf(float __complex__ z) {
float z_real = crealf(z);
float z_imag = cimagf(z);
if (__spirv_IsInf(z_real) && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLXF(__spirv_ocl_fabs(z_real), NAN);
return __SYCL_CMPLXF(__spirv_ocl_fabs(z_real),
std::numeric_limits<float>::quiet_NaN());
if (z_real == 0 && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLXF(NAN, z_real);
return __SYCL_CMPLXF(std::numeric_limits<float>::quiet_NaN(), z_real);
if (z_real == 0 && z_imag == 0)
return __SYCL_CMPLXF(1.0f, z_imag);
if (z_imag == 0 && !__spirv_IsFinite(z_real))
Expand All @@ -476,9 +496,10 @@ double __complex__ ccosh(double __complex__ z) {
double z_real = creal(z);
double z_imag = cimag(z);
if (__spirv_IsInf(z_real) && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLX(__spirv_ocl_fabs(z_real), NAN);
return __SYCL_CMPLX(__spirv_ocl_fabs(z_real),
std::numeric_limits<float>::quiet_NaN());
if (z_real == 0 && !__spirv_IsFinite(z_imag))
return __SYCL_CMPLX(NAN, z_real);
return __SYCL_CMPLX(std::numeric_limits<float>::quiet_NaN(), z_real);
if (z_real == 0 && z_imag == 0)
return __SYCL_CMPLX(1.0f, z_imag);
if (z_imag == 0 && !__spirv_IsFinite(z_real))
Expand Down Expand Up @@ -790,8 +811,9 @@ float __complex__ catanhf(float __complex__ z) {
return __SYCL_CMPLXF(__spirv_ocl_copysign(0.0f, z_real),
__spirv_ocl_copysign(__pi / 2.0f, z_imag));
if (__spirv_ocl_fabs(z_real) == 1.0f && z_imag == 0.0f)
return __SYCL_CMPLXF(__spirv_ocl_copysign(INFINITY, z_real),
__spirv_ocl_copysign(0.0f, z_imag));
return __SYCL_CMPLXF(
__spirv_ocl_copysign(std::numeric_limits<float>::infinity(), z_real),
__spirv_ocl_copysign(0.0f, z_imag));
float __complex__ t1 = 1.0f + z;
float __complex__ t2 = 1.0f - z;
float __complex__ t3 =
Expand Down Expand Up @@ -820,7 +842,7 @@ double __complex__ catanh(double __complex__ z) {
__spirv_ocl_copysign(__pi / 2.0, z_imag));
if (__spirv_ocl_fabs(z_real) == 1.0 && z_imag == 0.0)
return __SYCL_CMPLX(
__spirv_ocl_copysign(static_cast<double>(INFINITY), z_real),
__spirv_ocl_copysign(std::numeric_limits<double>::infinity(), z_real),
__spirv_ocl_copysign(0.0, z_imag));
double __complex__ t1 = 1.0 + z;
double __complex__ t2 = 1.0 - z;
Expand Down
Loading
Loading