Skip to content

[SYCL][NFC] Group all msvc specific math functions together #19913

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 2, 2025
Merged

[SYCL][NFC] Group all msvc specific math functions together #19913

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
820b447
[SYCL][NFC] Group all msvc specific math functions together
jinge90 Aug 28, 2025
02cd6bf
fix clang format
jinge90 Aug 28, 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
308 changes: 0 additions & 308 deletions libdevice/cmath_wrapper_fp64.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -181,312 +181,4 @@ double scalbn(double x, int exp) { return __devicelib_scalbn(x, exp); }

DEVICE_EXTERN_C_INLINE
double rint(double x) { return __spirv_ocl_rint(x); }

#if defined(_MSC_VER)
#include <math.h>
// FLOAT PROPERTIES
#define _D0 3 // little-endian, small long doubles
#define _D1 2
#define _D2 1
#define _D3 0

// IEEE 754 double properties
#define HUGE_EXP (int)(_DMAX * 900L / 1000)

#define NBITS (48 + _DOFF)

#define INIT(w0) \
{ 0, 0, 0, w0 }

// double declarations
union _Dval { // pun floating type as integer array
unsigned short _Sh[8];
double _Val;
};

union _Dconst { // pun float types as integer array
unsigned short _Word[4]; // TRANSITION, ABI: Twice as large as necessary.
double _Double;
};
#define DSIGN(x) (((_Dval *)(char *)&(x))->_Sh[_D0] & _DSIGN)

#define _Xbig (double)((NBITS + 1) * 347L / 1000)

DEVICE_EXTERN_C_INLINE
short _Dtest(double *px) { // categorize *px
_Dval *ps = (_Dval *)(char *)px;

short ret = 0;
if ((ps->_Sh[_D0] & _DMASK) == _DMAX << _DOFF) {
ret = (short)((ps->_Sh[_D0] & _DFRAC) != 0 || ps->_Sh[_D1] != 0 ||
ps->_Sh[_D2] != 0 || ps->_Sh[_D3] != 0
? _NANCODE
: _INFCODE);
} else if ((ps->_Sh[_D0] & ~_DSIGN) != 0 || ps->_Sh[_D1] != 0 ||
ps->_Sh[_D2] != 0 || ps->_Sh[_D3] != 0)
ret = (ps->_Sh[_D0] & _DMASK) == 0 ? _DENORM : _FINITE;

return ret;
}

// Returns _FP_LT, _FP_GT or _FP_EQ based on the ordering
// relationship between x and y.
DEVICE_EXTERN_C_INLINE
int _dpcomp(double x, double y) {
int res = 0;
if (_Dtest(&x) == _NANCODE || _Dtest(&y) == _NANCODE) {
// '0' means unordered.
return res;
}

if (x < y)
res |= _FP_LT;
else if (x > y)
res |= _FP_GT;
else
res |= _FP_EQ;

return res;
}

// Returns 0, if the sign bit is not set, and non-zero otherwise.
DEVICE_EXTERN_C_INLINE
int _dsign(double x) { return DSIGN(x); }

// fpclassify() equivalent with a pointer argument.
DEVICE_EXTERN_C_INLINE
short _dtest(double *px) {
switch (_Dtest(px)) {
case _DENORM:
return FP_SUBNORMAL;
case _FINITE:
return FP_NORMAL;
case _INFCODE:
return FP_INFINITE;
case _NANCODE:
return FP_NAN;
}

return FP_ZERO;
}

DEVICE_EXTERN_C_INLINE
short _Dnorm(_Dval *ps) { // normalize double fraction
short xchar;
unsigned short sign = (unsigned short)(ps->_Sh[_D0] & _DSIGN);

xchar = 1;
if ((ps->_Sh[_D0] &= _DFRAC) != 0 || ps->_Sh[_D1] || ps->_Sh[_D2] ||
ps->_Sh[_D3]) { // nonzero, scale
for (; ps->_Sh[_D0] == 0; xchar -= 16) { // shift left by 16
ps->_Sh[_D0] = ps->_Sh[_D1];
ps->_Sh[_D1] = ps->_Sh[_D2];
ps->_Sh[_D2] = ps->_Sh[_D3];
ps->_Sh[_D3] = 0;
}
for (; ps->_Sh[_D0] < 1 << _DOFF; --xchar) { // shift left by 1
ps->_Sh[_D0] = (unsigned short)(ps->_Sh[_D0] << 1 | ps->_Sh[_D1] >> 15);
ps->_Sh[_D1] = (unsigned short)(ps->_Sh[_D1] << 1 | ps->_Sh[_D2] >> 15);
ps->_Sh[_D2] = (unsigned short)(ps->_Sh[_D2] << 1 | ps->_Sh[_D3] >> 15);
ps->_Sh[_D3] <<= 1;
}
for (; 1 << (_DOFF + 1) <= ps->_Sh[_D0]; ++xchar) { // shift right by 1
ps->_Sh[_D3] = (unsigned short)(ps->_Sh[_D3] >> 1 | ps->_Sh[_D2] << 15);
ps->_Sh[_D2] = (unsigned short)(ps->_Sh[_D2] >> 1 | ps->_Sh[_D1] << 15);
ps->_Sh[_D1] = (unsigned short)(ps->_Sh[_D1] >> 1 | ps->_Sh[_D0] << 15);
ps->_Sh[_D0] >>= 1;
}
ps->_Sh[_D0] &= _DFRAC;
}
ps->_Sh[_D0] |= sign;
return xchar;
}

DEVICE_EXTERN_C_INLINE
short _Dscale(double *px, long lexp) { // scale *px by 2^xexp with checking
_Dval *ps = (_Dval *)(char *)px;
_Dconst _Inf = {INIT(_DMAX << _DOFF)};
short xchar = (short)((ps->_Sh[_D0] & _DMASK) >> _DOFF);

if (xchar == _DMAX)
return (short)((ps->_Sh[_D0] & _DFRAC) != 0 || ps->_Sh[_D1] != 0 ||
ps->_Sh[_D2] != 0 || ps->_Sh[_D3] != 0
? _NANCODE
: _INFCODE);
if (xchar == 0 && 0 < (xchar = _Dnorm(ps)))
return 0;

short ret = 0;
if (0 < lexp && _DMAX - xchar <= lexp) { // overflow, return +/-INF
*px = ps->_Sh[_D0] & _DSIGN ? -_Inf._Double : _Inf._Double;
ret = _INFCODE;
} else if (-xchar < lexp) { // finite result, repack
ps->_Sh[_D0] =
(unsigned short)(ps->_Sh[_D0] & ~_DMASK | (lexp + xchar) << _DOFF);
ret = _FINITE;
} else { // denormalized, scale
unsigned short sign = (unsigned short)(ps->_Sh[_D0] & _DSIGN);

ps->_Sh[_D0] = (unsigned short)(1 << _DOFF | ps->_Sh[_D0] & _DFRAC);
lexp += xchar - 1;
if (lexp < -(48 + 1 + _DOFF) ||
0 <= lexp) { // certain underflow, return +/-0
ps->_Sh[_D0] = sign;
ps->_Sh[_D1] = 0;
ps->_Sh[_D2] = 0;
ps->_Sh[_D3] = 0;
ret = 0;
} else { // nonzero, align fraction
short xexp = (short)lexp;
unsigned short psx = 0;
ret = _FINITE;

for (; xexp <= -16; xexp += 16) { // scale by words
psx = ps->_Sh[_D3] | (psx != 0 ? 1 : 0);
ps->_Sh[_D3] = ps->_Sh[_D2];
ps->_Sh[_D2] = ps->_Sh[_D1];
ps->_Sh[_D1] = ps->_Sh[_D0];
ps->_Sh[_D0] = 0;
}
if ((xexp = (short)-xexp) != 0) { // scale by bits
psx = (ps->_Sh[_D3] << (16 - xexp)) | (psx != 0 ? 1 : 0);
ps->_Sh[_D3] = (unsigned short)(ps->_Sh[_D3] >> xexp |
ps->_Sh[_D2] << (16 - xexp));
ps->_Sh[_D2] = (unsigned short)(ps->_Sh[_D2] >> xexp |
ps->_Sh[_D1] << (16 - xexp));
ps->_Sh[_D1] = (unsigned short)(ps->_Sh[_D1] >> xexp |
ps->_Sh[_D0] << (16 - xexp));
ps->_Sh[_D0] >>= xexp;
}

ps->_Sh[_D0] |= sign;
if ((0x8000 < psx || 0x8000 == psx && (ps->_Sh[_D3] & 0x0001) != 0) &&
(++ps->_Sh[_D3] & 0xffff) == 0 && (++ps->_Sh[_D2] & 0xffff) == 0 &&
(++ps->_Sh[_D1] & 0xffff) == 0)
++ps->_Sh[_D0]; // round up
else if (ps->_Sh[_D0] == sign && ps->_Sh[_D1] == 0 && ps->_Sh[_D2] == 0 &&
ps->_Sh[_D3] == 0)
ret = 0;
}
}
return ret;
}

DEVICE_EXTERN_C_INLINE
short _Exp(double *px, double y,
short eoff) { // compute y * e^(*px), (*px) finite, |y| not huge
static const double invln2 = 1.4426950408889634073599246810018921;
static const double c1 = 22713.0 / 32768.0;
static const double c2 = 1.4286068203094172321214581765680755e-6;
static const double p[] = {1.0, 420.30235984910635, 15132.70094680474802};
static const double q[] = {30.01511290683317, 3362.72154416553028,
30265.40189360949691};

_Dconst _Eps = {INIT((_DBIAS - NBITS - 1) << _DOFF)};
_Dconst _Inf = {INIT(_DMAX << _DOFF)};
short ret = 0;
if (*px < -HUGE_EXP || y == 0.0) // certain underflow
*px = __spirv_ocl_copysign(0.0, y);
else if (HUGE_EXP < *px) { // certain overflow
*px = __spirv_ocl_copysign(_Inf._Double, y);
ret = _INFCODE;
} else { // xexp won't overflow
double g = *px * invln2;
short xexp = (short)(g + (g < 0.0 ? -0.5 : +0.5));
g = xexp;
g = (*px - g * c1) - g * c2;
if (-_Eps._Double < g && g < _Eps._Double)
*px = y;
else { // g * g worth computing
const double z = g * g;
const double w = (q[0] * z + q[1]) * z + q[2];

g *= (z + p[1]) * z + p[2];
*px = (w + g) / (w - g) * 2.0 * y;
--xexp;
}
ret = _Dscale(px, (long)xexp + eoff);
}

return ret;
}

DEVICE_EXTERN_C_INLINE
double _Cosh(double x, double y) { // compute y * cosh(x), |y| <= 1
switch (_Dtest(&x)) { // test for special codes
case _NANCODE:
case _INFCODE:
return x;
case 0:
return y;
default: // finite
if (y == 0.0)
return y;

if (x < 0.0)
x = -x;

if (x < _Xbig) { // worth adding in exp(-x)
_Exp(&x, 1.0, -1);
return y * (x + 0.25 / x);
}
_Exp(&x, y, -1);
return x;
}
}

DEVICE_EXTERN_C_INLINE
double _Poly(double x, const double *tab, int n) { // compute polynomial
double y;

for (y = *tab; 0 <= --n;)
y = y * x + *++tab;

return y;
}

DEVICE_EXTERN_C_INLINE
double _Sinh(double x, double y) { // compute y * sinh(x), |y| <= 1

short neg;
// coefficients
static const double p[] = {0.0000000001632881, 0.0000000250483893,
0.0000027557344615, 0.0001984126975233,
0.0083333333334816, 0.1666666666666574,
1.0000000000000001};
_Dconst _Rteps = {INIT((_DBIAS - NBITS / 2) << _DOFF)};
switch (_Dtest(&x)) { // test for special codes
case _NANCODE:
return x;
case _INFCODE:
return y != 0.0 ? x : DSIGN(x) ? -y : y;
case 0:
return x * y;
default: // finite
if (y == 0.0)
return x < 0.0 ? -y : y;

if (x < 0.0) {
x = -x;
neg = 1;
} else
neg = 0;

if (x < _Rteps._Double)
x *= y; // x tiny
else if (x < 1.0) {
double w = x * x;

x += x * w * _Poly(w, p, 5);
x *= y;
} else if (x < _Xbig) { // worth adding in exp(-x)
_Exp(&x, 1.0, -1);
x = y * (x - 0.25 / x);
} else
_Exp(&x, y, -1);

return neg ? -x : x;
}
}
#endif // defined(_WIN32)
#endif // __SPIR__ || __SPIRV__
Loading