Alexandre Julliard : msvcrt: Use the jn()/yn() implementation from the bundled musl library.
Module: wine Branch: master Commit: ec31d30ec09efe29c5dbce6ab044c0989b9fc153 URL: https://gitlab.winehq.org/wine/wine/-/commit/ec31d30ec09efe29c5dbce6ab044c09...
Author: Alexandre Julliard julliard@winehq.org Date: Tue Apr 4 16:31:20 2023 +0200
msvcrt: Use the jn()/yn() implementation from the bundled musl library.
---
dlls/msvcrt/math.c | 200 ------------------------------------------------ libs/musl/src/math/jn.c | 2 +- 2 files changed, 1 insertion(+), 201 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 4437cbbf48e..5ff75274dbc 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -5510,206 +5510,6 @@ int CDECL _isnan(double num) return (u.i & ~0ull >> 1) > 0x7ffull << 52; }
-/********************************************************************* - * _jn (MSVCRT.@) - * - * Copied from musl: src/math/jn.c - */ -double CDECL _jn(int n, double x) -{ - static const double invsqrtpi = 5.64189583547756279280e-01; - - unsigned int ix, lx; - int nm1, i, sign; - double a, b, temp; - - ix = *(ULONGLONG*)&x >> 32; - lx = *(ULONGLONG*)&x; - sign = ix >> 31; - ix &= 0x7fffffff; - - if ((ix | (lx | -lx) >> 31) > 0x7ff00000) /* nan */ - return x; - - if (n == 0) - return _j0(x); - if (n < 0) { - nm1 = -(n + 1); - x = -x; - sign ^= 1; - } else { - nm1 = n-1; - } - if (nm1 == 0) - return j1(x); - - sign &= n; /* even n: 0, odd n: signbit(x) */ - x = fabs(x); - if ((ix | lx) == 0 || ix == 0x7ff00000) /* if x is 0 or inf */ - b = 0.0; - else if (nm1 < x) { - if (ix >= 0x52d00000) { /* x > 2**302 */ - switch(nm1 & 3) { - case 0: - temp = -cos(x) + sin(x); - break; - case 1: - temp = -cos(x) - sin(x); - break; - case 2: - temp = cos(x) - sin(x); - break; - default: - temp = cos(x) + sin(x); - break; - } - b = invsqrtpi * temp / sqrt(x); - } else { - a = _j0(x); - b = _j1(x); - for (i = 0; i < nm1; ) { - i++; - temp = b; - b = b * (2.0 * i / x) - a; /* avoid underflow */ - a = temp; - } - } - } else { - if (ix < 0x3e100000) { /* x < 2**-29 */ - if (nm1 > 32) /* underflow */ - b = 0.0; - else { - temp = x * 0.5; - b = temp; - a = 1.0; - for (i = 2; i <= nm1 + 1; i++) { - a *= (double)i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b / a; - } - } else { - double t, q0, q1, w, h, z, tmp, nf; - int k; - - nf = nm1 + 1.0; - w = 2 * nf / x; - h = 2 / x; - z = w + h; - q0 = w; - q1 = w * z - 1.0; - k = 1; - while (q1 < 1.0e9) { - k += 1; - z += h; - tmp = z * q1 - q0; - q0 = q1; - q1 = tmp; - } - for (t = 0.0, i = k; i >= 0; i--) - t = 1 / (2 * (i + nf) / x - t); - a = t; - b = 1.0; - tmp = nf * log(fabs(w)); - if (tmp < 7.09782712893383973096e+02) { - for (i = nm1; i > 0; i--) { - temp = b; - b = b * (2.0 * i) / x - a; - a = temp; - } - } else { - for (i = nm1; i > 0; i--) { - temp = b; - b = b * (2.0 * i) / x - a; - a = temp; - /* scale b to avoid spurious overflow */ - if (b > 0x1p500) { - a /= b; - t /= b; - b = 1.0; - } - } - } - z = j0(x); - w = j1(x); - if (fabs(z) >= fabs(w)) - b = t * z / b; - else - b = t * w / a; - } - } - return sign ? -b : b; -} - -/********************************************************************* - * _yn (MSVCRT.@) - * - * Copied from musl: src/math/jn.c - */ -double CDECL _yn(int n, double x) -{ - static const double invsqrtpi = 5.64189583547756279280e-01; - - unsigned int ix, lx, ib; - int nm1, sign, i; - double a, b, temp; - - ix = *(ULONGLONG*)&x >> 32; - lx = *(ULONGLONG*)&x; - sign = ix >> 31; - ix &= 0x7fffffff; - - if ((ix | (lx | -lx) >> 31) > 0x7ff00000) /* nan */ - return x; - if (sign && (ix | lx) != 0) /* x < 0 */ - return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x)); - if (ix == 0x7ff00000) - return 0.0; - - if (n == 0) - return y0(x); - if (n < 0) { - nm1 = -(n + 1); - sign = n & 1; - } else { - nm1 = n - 1; - sign = 0; - } - if (nm1 == 0) - return sign ? -y1(x) : y1(x); - - if (ix >= 0x52d00000) { /* x > 2**302 */ - switch(nm1 & 3) { - case 0: - temp = -sin(x) - cos(x); - break; - case 1: - temp = -sin(x) + cos(x); - break; - case 2: - temp = sin(x) + cos(x); - break; - default: - temp = sin(x) - cos(x); - break; - } - b = invsqrtpi * temp / sqrt(x); - } else { - a = y0(x); - b = y1(x); - /* quit if b is -inf */ - ib = *(ULONGLONG*)&b >> 32; - for (i = 0; i < nm1 && ib != 0xfff00000;) { - i++; - temp = b; - b = (2.0 * i / x) * b - a; - ib = *(ULONGLONG*)&b >> 32; - a = temp; - } - } - return sign ? -b : b; -} - #if _MSVCR_VER>=120
/********************************************************************* diff --git a/libs/musl/src/math/jn.c b/libs/musl/src/math/jn.c index 55c05157be0..f7e679e25ea 100644 --- a/libs/musl/src/math/jn.c +++ b/libs/musl/src/math/jn.c @@ -225,7 +225,7 @@ double __cdecl _yn(int n, double x) if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */ return x; if (sign && (ix|lx)!=0) /* x < 0 */ - return 0/0.0; + return math_error(_DOMAIN, "_y1", x, 0, 0 / (x - x)); if (ix == 0x7ff00000) return 0.0;
-
Alexandre Julliard