Skip to content

Commit 781898e

Browse files
oleksandr-pavlykoleksandr-pavlyk
committed
Untangle special case handling, hoping to fix SYCL OS build
1 parent 0b4fc0c commit 781898e

File tree

10 files changed

+214
-156
lines changed

10 files changed

+214
-156
lines changed

‎dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp

Lines changed: 25 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -67,47 +67,52 @@ template <typename argT, typename resT> struct AcosFunctor
6767
{
6868
if constexpr (is_complex<argT>::value) {
6969
using realT = typename argT::value_type;
70+
71+
constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();
72+
7073
const realT x = std::real(in);
7174
const realT y = std::imag(in);
7275

73-
if (std::isnan(x) || std::isnan(y)) {
74-
/* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
75-
if (std::isinf(x)) {
76-
return resT{y + y, -std::numeric_limits<realT>::infinity()};
77-
}
76+
if (std::isnan(x)) {
7877
/* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
7978
if (std::isinf(y)) {
80-
return resT{x + x, -y};
79+
return resT{q_nan, -y};
80+
}
81+
82+
/* all other cases involving NaN return NaN + I*NaN. */
83+
return resT{q_nan, q_nan};
84+
}
85+
if (std::isnan(y)) {
86+
/* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
87+
if (std::isinf(x)) {
88+
return resT{q_nan, -std::numeric_limits<realT>::infinity()};
8189
}
8290
/* acos(0 + I*NaN) = PI/2 + I*NaN with inexact */
8391
if (x == realT(0)) {
8492
const realT res_re = std::atan(realT(1)) * 2; // PI/2
85-
return resT{res_re, y + y};
93+
return resT{res_re, q_nan};
8694
}
87-
/*
88-
* All other cases involving NaN return NaN + I*NaN.
89-
* C99 leaves it optional whether to raise invalid if one of
90-
* the arguments is not NaN, so we opt not to raise it.
91-
*/
92-
return resT{std::numeric_limits<realT>::quiet_NaN(),
93-
std::numeric_limits<realT>::quiet_NaN()};
95+
96+
/* all other cases involving NaN return NaN + I*NaN. */
97+
return resT{q_nan, q_nan};
9498
}
99+
95100
/*
96101
* For large x or y including acos(+-Inf + I*+-Inf)
97102
*/
98-
const realT RECIP_EPSILON =
103+
constexpr realT r_eps =
99104
realT(1) / std::numeric_limits<realT>::epsilon();
100-
if (std::abs(x) > RECIP_EPSILON || std::abs(y) > RECIP_EPSILON) {
105+
if (std::abs(x) > r_eps || std::abs(y) > r_eps) {
101106
argT log_in = std::log(in);
107+
102108
const realT wx = std::real(log_in);
103109
const realT wy = std::imag(log_in);
104110
const realT rx = std::abs(wy);
111+
105112
realT ry = wx + std::log(realT(2));
106-
if (!std::signbit(y)) {
107-
ry = -ry;
108-
}
109-
return resT{rx, ry};
113+
return resT{rx, (std::signbit(y)) ? ry : -ry};
110114
}
115+
111116
/* ordinary cases */
112117
return std::acos(in);
113118
}

‎dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp

Lines changed: 26 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -67,54 +67,54 @@ template <typename argT, typename resT> struct AcoshFunctor
6767
{
6868
if constexpr (is_complex<argT>::value) {
6969
using realT = typename argT::value_type;
70+
71+
constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();
7072
/*
7173
* acosh(in) = I*acos(in) or -I*acos(in)
7274
* where the sign is chosen so Re(acosh(in)) >= 0.
7375
* So, we first calculate acos(in) and then acosh(in).
7476
*/
7577
const realT x = std::real(in);
7678
const realT y = std::imag(in);
77-
const realT RECIP_EPSILON =
78-
realT(1) / std::numeric_limits<realT>::epsilon();
79+
7980
resT acos_in;
80-
if (std::isnan(x) || std::isnan(y)) {
81+
if (std::isnan(x)) {
82+
/* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
83+
if (std::isinf(y)) {
84+
acos_in = resT{q_nan, -y};
85+
}
86+
else {
87+
acos_in = resT{q_nan, q_nan};
88+
}
89+
}
90+
else if (std::isnan(y)) {
8191
/* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
92+
constexpr realT inf = std::numeric_limits<realT>::infinity();
93+
8294
if (std::isinf(x)) {
83-
acos_in =
84-
resT{y + y, -std::numeric_limits<realT>::infinity()};
85-
}
86-
/* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
87-
else if (std::isinf(y)) {
88-
acos_in = resT{x + x, -y};
95+
acos_in = resT{q_nan, -inf};
8996
}
90-
/* acos(0 + I*NaN) = PI/2 + I*NaN with inexact */
97+
/* acos(0 + I*NaN) = Pi/2 + I*NaN with inexact */
9198
else if (x == realT(0)) {
92-
const realT res_re = std::atan(1) * 2; // PI/2
93-
acos_in = resT{res_re, y + y};
99+
const realT pi_half = std::atan(1) * 2;
100+
acos_in = resT{pi_half, q_nan};
94101
}
95-
/*
96-
* All other cases involving NaN return NaN + I*NaN.
97-
* C99 leaves it optional whether to raise invalid if one of
98-
* the arguments is not NaN, so we opt not to raise it.
99-
*/
100102
else {
101-
acos_in = resT{std::numeric_limits<realT>::quiet_NaN(),
102-
std::numeric_limits<realT>::quiet_NaN()};
103+
acos_in = resT{q_nan, q_nan};
103104
}
104105
}
106+
107+
constexpr realT r_eps =
108+
realT(1) / std::numeric_limits<realT>::epsilon();
105109
/*
106110
* For large x or y including acos(+-Inf + I*+-Inf)
107111
*/
108-
else if (std::abs(x) > RECIP_EPSILON || std::abs(y) > RECIP_EPSILON)
109-
{
112+
if (std::abs(x) > r_eps || std::abs(y) > r_eps) {
110113
const realT wx = std::real(std::log(in));
111114
const realT wy = std::imag(std::log(in));
112115
const realT rx = std::abs(wy);
113116
realT ry = wx + std::log(realT(2));
114-
if (!std::signbit(y)) {
115-
ry = -ry;
116-
}
117-
acos_in = resT{rx, ry};
117+
acos_in = resT{rx, (std::signbit(y)) ? ry : -ry};
118118
}
119119
else {
120120
/* ordinary cases */
@@ -124,6 +124,7 @@ template <typename argT, typename resT> struct AcoshFunctor
124124
/* Now we calculate acosh(z) */
125125
const realT rx = std::real(acos_in);
126126
const realT ry = std::imag(acos_in);
127+
127128
/* acosh(NaN + I*NaN) = NaN + I*NaN */
128129
if (std::isnan(rx) && std::isnan(ry)) {
129130
return resT{ry, rx};

‎dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp

Lines changed: 29 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,9 @@ template <typename argT, typename resT> struct AsinFunctor
6767
{
6868
if constexpr (is_complex<argT>::value) {
6969
using realT = typename argT::value_type;
70+
71+
constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();
72+
7073
/*
7174
* asin(in) = I * conj( asinh(I * conj(in)) )
7275
* so we first calculate w = asinh(I * conj(in)) with
@@ -77,63 +80,61 @@ template <typename argT, typename resT> struct AsinFunctor
7780
const realT x = std::imag(in);
7881
const realT y = std::real(in);
7982

80-
if (std::isnan(x) || std::isnan(y)) {
81-
/* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */
82-
if (std::isinf(x)) {
83-
const realT asinh_re = x;
84-
const realT asinh_im = y + y;
85-
return resT{asinh_im, asinh_re};
86-
}
83+
if (std::isnan(x)) {
8784
/* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
8885
if (std::isinf(y)) {
8986
const realT asinh_re = y;
90-
const realT asinh_im = x + x;
87+
const realT asinh_im = q_nan;
9188
return resT{asinh_im, asinh_re};
9289
}
9390
/* asinh(NaN + I*0) = NaN + I*0 */
9491
if (y == realT(0)) {
95-
const realT asinh_re = x + x;
92+
const realT asinh_re = q_nan;
9693
const realT asinh_im = y;
9794
return resT{asinh_im, asinh_re};
9895
}
99-
/*
100-
* All other cases involving NaN return NaN + I*NaN.
101-
* C99 leaves it optional whether to raise invalid if one of
102-
* the arguments is not NaN, so we opt not to raise it.
103-
*/
104-
return resT{std::numeric_limits<realT>::quiet_NaN(),
105-
std::numeric_limits<realT>::quiet_NaN()};
96+
/* All other cases involving NaN return NaN + I*NaN. */
97+
return resT{q_nan, q_nan};
10698
}
99+
else if (std::isnan(y)) {
100+
/* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */
101+
if (std::isinf(x)) {
102+
const realT asinh_re = x;
103+
const realT asinh_im = q_nan;
104+
return resT{asinh_im, asinh_re};
105+
}
106+
/* All other cases involving NaN return NaN + I*NaN. */
107+
return resT{q_nan, q_nan};
108+
}
109+
107110
/*
108111
* For large x or y including asinh(+-Inf + I*+-Inf)
109112
* asinh(in) = sign(x)*log(sign(x)*in) + O(1/in^2) as in ->
110113
* infinity The above formula works for the imaginary part as well,
111114
* because Im(asinh(in)) = sign(x)*atan2(sign(x)*y, fabs(x)) +
112115
* O(y/in^3) as in -> infinity, uniformly in y
113116
*/
114-
const realT RECIP_EPSILON =
117+
constexpr realT r_eps =
115118
realT(1) / std::numeric_limits<realT>::epsilon();
116-
const resT z = {x, y};
117-
if (std::abs(x) > RECIP_EPSILON || std::abs(y) > RECIP_EPSILON) {
119+
if (std::abs(x) > r_eps || std::abs(y) > r_eps) {
120+
const resT z = {x, y};
118121
realT wx, wy;
119122
if (!std::signbit(x)) {
120-
wx = std::real(std::log(z));
121-
wy = std::imag(std::log(z));
122-
wx += std::log(realT(2));
123+
auto log_z = std::log(z);
124+
wx = std::real(log_z) + std::log(realT(2));
125+
wy = std::imag(log_z);
123126
}
124127
else {
125-
wx = std::real(std::log(-z));
126-
wy = std::imag(std::log(-z));
127-
wx += std::log(realT(2));
128+
auto log_mz = std::log(-z);
129+
wx = std::real(log_mz) + std::log(realT(2));
130+
wy = std::imag(log_mz);
128131
}
129132
const realT asinh_re = std::copysign(wx, x);
130133
const realT asinh_im = std::copysign(wy, y);
131134
return resT{asinh_im, asinh_re};
132135
}
133136
/* ordinary cases */
134-
const realT asinh_re = std::real(std::asinh(z));
135-
const realT asinh_im = std::imag(std::asinh(z));
136-
return resT{asinh_im, asinh_re};
137+
return std::asin(in);
137138
}
138139
else {
139140
static_assert(std::is_floating_point_v<argT> ||

‎dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp

Lines changed: 22 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -67,28 +67,32 @@ template <typename argT, typename resT> struct AsinhFunctor
6767
{
6868
if constexpr (is_complex<argT>::value) {
6969
using realT = typename argT::value_type;
70+
71+
constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();
72+
7073
const realT x = std::real(in);
7174
const realT y = std::imag(in);
72-
if (std::isnan(x) || std::isnan(y)) {
73-
/* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */
74-
if (std::isinf(x)) {
75-
return resT{x, y + y};
76-
}
75+
76+
if (std::isnan(x)) {
7777
/* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
7878
if (std::isinf(y)) {
79-
return resT{y, x + x};
79+
return resT{y, q_nan};
8080
}
8181
/* asinh(NaN + I*0) = NaN + I*0 */
8282
if (y == realT(0)) {
83-
return resT{x + x, y};
83+
return resT{q_nan, y};
8484
}
85-
/*
86-
* All other cases involving NaN return NaN + I*NaN.
87-
* C99 leaves it optional whether to raise invalid if one of
88-
* the arguments is not NaN, so we opt not to raise it.
89-
*/
90-
return resT{std::numeric_limits<realT>::quiet_NaN(),
91-
std::numeric_limits<realT>::quiet_NaN()};
85+
/* All other cases involving NaN return NaN + I*NaN. */
86+
return resT{q_nan, q_nan};
87+
}
88+
89+
if (std::isnan(y)) {
90+
/* asinh(+-Inf + I*NaN) = +-Inf + I*NaN */
91+
if (std::isinf(x)) {
92+
return resT{x, q_nan};
93+
}
94+
/* All other cases involving NaN return NaN + I*NaN. */
95+
return resT{q_nan, q_nan};
9296
}
9397

9498
/*
@@ -98,16 +102,18 @@ template <typename argT, typename resT> struct AsinhFunctor
98102
* because Im(asinh(in)) = sign(x)*atan2(sign(x)*y, fabs(x)) +
99103
* O(y/in^3) as in -> infinity, uniformly in y
100104
*/
101-
const realT RECIP_EPSILON =
105+
constexpr realT r_eps =
102106
realT(1) / std::numeric_limits<realT>::epsilon();
103-
if (std::abs(x) > RECIP_EPSILON || std::abs(y) > RECIP_EPSILON) {
107+
108+
if (std::abs(x) > r_eps || std::abs(y) > r_eps) {
104109
resT log_in = (std::signbit(x)) ? std::log(-in) : std::log(in);
105110
realT wx = std::real(log_in) + std::log(realT(2));
106111
realT wy = std::imag(log_in);
107112
const realT res_re = std::copysign(wx, x);
108113
const realT res_im = std::copysign(wy, y);
109114
return resT{res_re, res_im};
110115
}
116+
111117
/* ordinary cases */
112118
return std::asinh(in);
113119
}

0 commit comments

Comments
 (0)