diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 9deaf832d15..e943076c1c7 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -311,6 +311,244 @@ float CDECL _logbf(float x) #endif +/* Copied from musl: src/math/scalbn.c */ +static double __scalbn(double x, int n) +{ + union {double f; UINT64 i;} u; + double y = x; + + if (n > 1023) { + y *= 0x1p1023; + n -= 1023; + if (n > 1023) { + y *= 0x1p1023; + n -= 1023; + if (n > 1023) + n = 1023; + } + } else if (n < -1022) { + /* make sure final n < -53 to avoid double + rounding in the subnormal range */ + y *= 0x1p-1022 * 0x1p53; + n += 1022 - 53; + if (n < -1022) { + y *= 0x1p-1022 * 0x1p53; + n += 1022 - 53; + if (n < -1022) + n = -1022; + } + } + u.i = (UINT64)(0x3ff + n) << 52; + x = y * u.f; + return x; +} + +/* Copied from musl: src/math/__rem_pio2_large.c */ +static int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec) +{ + static const int init_jk[] = {3, 4}; + static const INT32 ipio2[] = { + 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, + 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, + 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, + 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, + 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, + 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, + 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, + 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, + 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, + 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, + 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, + }; + static const double PIo2[] = { + 1.57079625129699707031e+00, + 7.54978941586159635335e-08, + 5.39030252995776476554e-15, + 3.28200341580791294123e-22, + 1.27065575308067607349e-29, + 1.22933308981111328932e-36, + 2.73370053816464559624e-44, + 2.16741683877804819444e-51, + }; + + INT32 jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; + double z, fw, f[20], fq[20] = {0}, q[20]; + + /* initialize jk*/ + jk = init_jk[prec]; + jp = jk; + + /* determine jx,jv,q0, note that 3>q0 */ + jx = nx - 1; + jv = (e0 - 3) / 24; + if(jv < 0) jv = 0; + q0 = e0 - 24 * (jv + 1); + + /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ + j = jv - jx; + m = jx + jk; + for (i = 0; i <= m; i++, j++) + f[i] = j < 0 ? 0.0 : (double)ipio2[j]; + + /* compute q[0],q[1],...q[jk] */ + for (i = 0; i <= jk; i++) { + for (j = 0, fw = 0.0; j <= jx; j++) + fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + + jz = jk; +recompute: + /* distill q[] into iq[] reversingly */ + for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) { + fw = (double)(INT32)(0x1p-24 * z); + iq[i] = (INT32)(z - 0x1p24 * fw); + z = q[j - 1] + fw; + } + + /* compute n */ + z = __scalbn(z, q0); /* actual value of z */ + z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */ + n = (INT32)z; + z -= (double)n; + ih = 0; + if (q0 > 0) { /* need iq[jz-1] to determine n */ + i = iq[jz - 1] >> (24 - q0); + n += i; + iq[jz - 1] -= i << (24 - q0); + ih = iq[jz - 1] >> (23 - q0); + } + else if (q0 == 0) ih = iq[jz - 1] >> 23; + else if (z >= 0.5) ih = 2; + + if (ih > 0) { /* q > 0.5 */ + n += 1; + carry = 0; + for (i = 0; i < jz; i++) { /* compute 1-q */ + j = iq[i]; + if (carry == 0) { + if (j != 0) { + carry = 1; + iq[i] = 0x1000000 - j; + } + } else + iq[i] = 0xffffff - j; + } + if (q0 > 0) { /* rare case: chance is 1 in 12 */ + switch(q0) { + case 1: + iq[jz - 1] &= 0x7fffff; + break; + case 2: + iq[jz - 1] &= 0x3fffff; + break; + } + } + if (ih == 2) { + z = 1.0 - z; + if (carry != 0) + z -= __scalbn(1.0, q0); + } + } + + /* check if recomputation is needed */ + if (z == 0.0) { + j = 0; + for (i = jz - 1; i >= jk; i--) j |= iq[i]; + if (j == 0) { /* need recomputation */ + for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */ + + for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */ + f[jx + i] = (double)ipio2[jv + i]; + for (j = 0, fw = 0.0; j <= jx; j++) + fw += x[j] * f[jx + i - j]; + q[i] = fw; + } + jz += k; + goto recompute; + } + } + + /* chop off zero terms */ + if (z == 0.0) { + jz -= 1; + q0 -= 24; + while (iq[jz] == 0) { + jz--; + q0 -= 24; + } + } else { /* break z into 24-bit if necessary */ + z = __scalbn(z, -q0); + if (z >= 0x1p24) { + fw = (double)(INT32)(0x1p-24 * z); + iq[jz] = (INT32)(z - 0x1p24 * fw); + jz += 1; + q0 += 24; + iq[jz] = (INT32)fw; + } else + iq[jz] = (INT32)z; + } + + /* convert integer "bit" chunk to floating-point value */ + fw = __scalbn(1.0, q0); + for (i = jz; i >= 0; i--) { + q[i] = fw * (double)iq[i]; + fw *= 0x1p-24; + } + + /* compute PIo2[0,...,jp]*q[jz,...,0] */ + for(i = jz; i >= 0; i--) { + for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) + fw += PIo2[k] * q[i + k]; + fq[jz - i] = fw; + } + + /* compress fq[] into y[] */ + switch(prec) { + case 0: + fw = 0.0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + y[0] = ih == 0 ? fw : -fw; + break; + case 1: + case 2: + fw = 0.0; + for (i = jz; i >= 0; i--) + fw += fq[i]; + fw = (double)fw; + y[0] = ih==0 ? fw : -fw; + fw = fq[0] - fw; + for (i = 1; i <= jz; i++) + fw += fq[i]; + y[1] = ih == 0 ? fw : -fw; + break; + case 3: /* painful */ + for (i = jz; i > 0; i--) { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (i = jz; i > 1; i--) { + fw = fq[i - 1] + fq[i]; + fq[i] += fq[i - 1] - fw; + fq[i - 1] = fw; + } + for (fw = 0.0, i = jz; i >= 2; i--) + fw += fq[i]; + if (ih == 0) { + y[0] = fq[0]; + y[1] = fq[1]; + y[2] = fw; + } else { + y[0] = -fq[0]; + y[1] = -fq[1]; + y[2] = -fw; + } + } + return n & 7; +} + #ifndef __i386__ /********************************************************************* @@ -609,14 +847,151 @@ float CDECL atan2f( float y, float x ) } } +/* Copied from musl: src/math/__rem_pio2f.c */ +static int __rem_pio2f(float x, double *y) +{ + static const double toint = 1.5 / DBL_EPSILON, + pio4 = 0x1.921fb6p-1, + invpio2 = 6.36619772367581382433e-01, + pio2_1 = 1.57079631090164184570e+00, + pio2_1t = 1.58932547735281966916e-08; + + union {float f; uint32_t i;} u = {x}; + double tx[1], ty[1], fn; + UINT32 ix; + int n, sign, e0; + + ix = u.i & 0x7fffffff; + /* 25+53 bit pi is good enough for medium size */ + if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ + /* Use a specialized rint() to get fn. */ + fn = fp_barrier(x * invpio2 + toint) - toint; + n = (int)fn; + *y = x - fn * pio2_1 - fn * pio2_1t; + /* Matters with directed rounding. */ + if (*y < -pio4) { + n--; + fn--; + *y = x - fn * pio2_1 - fn * pio2_1t; + } else if (*y > pio4) { + n++; + fn++; + *y = x - fn * pio2_1 - fn * pio2_1t; + } + return n; + } + if(ix >= 0x7f800000) { /* x is inf or NaN */ + *y = x - x; + return 0; + } + /* scale x into [2^23, 2^24-1] */ + sign = u.i >> 31; + e0 = (ix >> 23) - (0x7f + 23); /* e0 = ilogb(|x|)-23, positive */ + u.i = ix - (e0 << 23); + tx[0] = u.f; + n = __rem_pio2_large(tx, ty, e0, 1, 0); + if (sign) { + *y = -ty[0]; + return -n; + } + *y = ty[0]; + return n; +} + +/* Copied from musl: src/math/__sindf.c */ +static float __sindf(double x) +{ + static const double S1 = -0x15555554cbac77.0p-55, + S2 = 0x111110896efbb2.0p-59, + S3 = -0x1a00f9e2cae774.0p-65, + S4 = 0x16cd878c3b46a7.0p-71; + + double r, s, w, z; + + z = x * x; + w = z * z; + r = S3 + z * S4; + s = z * x; + return (x + s * (S1 + z * S2)) + s * w * r; +} + +/* Copied from musl: src/math/__cosdf.c */ +static float __cosdf(double x) +{ + static const double C0 = -0x1ffffffd0c5e81.0p-54, + C1 = 0x155553e1053a42.0p-57, + C2 = -0x16c087e80f1e27.0p-62, + C3 = 0x199342e0ee5069.0p-68; + double r, w, z; + + z = x * x; + w = z * z; + r = C2 + z * C3; + return ((1.0 + z * C0) + w * C1) + (w * z) * r; +} + /********************************************************************* * cosf (MSVCRT.@) + * + * Copied from musl: src/math/cosf.c */ float CDECL cosf( float x ) { - float ret = unix_funcs->cosf( x ); - if (!isfinite(x)) return math_error(_DOMAIN, "cosf", x, 0, ret); - return ret; + static const double c1pio2 = 1*M_PI_2, + c2pio2 = 2*M_PI_2, + c3pio2 = 3*M_PI_2, + c4pio2 = 4*M_PI_2; + + double y; + UINT32 ix; + unsigned n, sign; + + ix = *(UINT32*)&x; + sign = ix >> 31; + ix &= 0x7fffffff; + + if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + /* raise inexact if x != 0 */ + fp_barrierf(x + 0x1p120f); + return 1.0f; + } + return __cosdf(x); + } + if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */ + if (ix > 0x4016cbe3) /* |x| ~> 3*pi/4 */ + return -__cosdf(sign ? x + c2pio2 : x - c2pio2); + else { + if (sign) + return __sindf(x + c1pio2); + else + return __sindf(c1pio2 - x); + } + } + if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */ + if (ix > 0x40afeddf) /* |x| ~> 7*pi/4 */ + return __cosdf(sign ? x + c4pio2 : x - c4pio2); + else { + if (sign) + return __sindf(-x - c3pio2); + else + return __sindf(x - c3pio2); + } + } + + /* cos(Inf or NaN) is NaN */ + if (isinf(x)) return math_error(_DOMAIN, "cosf", x, 0, x - x); + if (ix >= 0x7f800000) + return x - x; + + /* general argument reduction needed */ + n = __rem_pio2f(x, &y); + switch (n & 3) { + case 0: return __cosdf(y); + case 1: return __sindf(-y); + case 2: return -__cosdf(y); + default: return __sindf(y); + } } /********************************************************************* @@ -1343,38 +1718,6 @@ double CDECL atan2( double y, double x ) } } -/* Copied from musl: src/math/scalbn.c */ -static double __scalbn(double x, int n) -{ - union {double f; UINT64 i;} u; - double y = x; - - if (n > 1023) { - y *= 0x1p1023; - n -= 1023; - if (n > 1023) { - y *= 0x1p1023; - n -= 1023; - if (n > 1023) - n = 1023; - } - } else if (n < -1022) { - /* make sure final n < -53 to avoid double - rounding in the subnormal range */ - y *= 0x1p-1022 * 0x1p53; - n += 1022 - 53; - if (n < -1022) { - y *= 0x1p-1022 * 0x1p53; - n += 1022 - 53; - if (n < -1022) - n = -1022; - } - } - u.i = (UINT64)(0x3ff + n) << 52; - x = y * u.f; - return x; -} - /* Copied from musl: src/math/rint.c */ static double __rint(double x) { @@ -1402,212 +1745,6 @@ static double __rint(double x) return y; } -/* Copied from musl: src/math/__rem_pio2_large.c */ -static int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec) -{ - static const int init_jk[] = {3, 4}; - static const INT32 ipio2[] = { - 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, - 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, - 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, - 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, - 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, - 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, - 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, - 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, - 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, - 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, - 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, - }; - static const double PIo2[] = { - 1.57079625129699707031e+00, - 7.54978941586159635335e-08, - 5.39030252995776476554e-15, - 3.28200341580791294123e-22, - 1.27065575308067607349e-29, - 1.22933308981111328932e-36, - 2.73370053816464559624e-44, - 2.16741683877804819444e-51, - }; - - INT32 jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; - double z, fw, f[20], fq[20], q[20]; - - /* initialize jk*/ - jk = init_jk[prec]; - jp = jk; - - /* determine jx,jv,q0, note that 3>q0 */ - jx = nx - 1; - jv = (e0 - 3) / 24; - if(jv < 0) jv = 0; - q0 = e0 - 24 * (jv + 1); - - /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ - j = jv - jx; - m = jx + jk; - for (i = 0; i <= m; i++, j++) - f[i] = j < 0 ? 0.0 : (double)ipio2[j]; - - /* compute q[0],q[1],...q[jk] */ - for (i = 0; i <= jk; i++) { - for (j = 0, fw = 0.0; j <= jx; j++) - fw += x[j] * f[jx + i - j]; - q[i] = fw; - } - - jz = jk; -recompute: - /* distill q[] into iq[] reversingly */ - for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) { - fw = (double)(INT32)(0x1p-24 * z); - iq[i] = (INT32)(z - 0x1p24 * fw); - z = q[j - 1] + fw; - } - - /* compute n */ - z = __scalbn(z, q0); /* actual value of z */ - z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */ - n = (INT32)z; - z -= (double)n; - ih = 0; - if (q0 > 0) { /* need iq[jz-1] to determine n */ - i = iq[jz - 1] >> (24 - q0); - n += i; - iq[jz - 1] -= i << (24 - q0); - ih = iq[jz - 1] >> (23 - q0); - } - else if (q0 == 0) ih = iq[jz - 1] >> 23; - else if (z >= 0.5) ih = 2; - - if (ih > 0) { /* q > 0.5 */ - n += 1; - carry = 0; - for (i = 0; i < jz; i++) { /* compute 1-q */ - j = iq[i]; - if (carry == 0) { - if (j != 0) { - carry = 1; - iq[i] = 0x1000000 - j; - } - } else - iq[i] = 0xffffff - j; - } - if (q0 > 0) { /* rare case: chance is 1 in 12 */ - switch(q0) { - case 1: - iq[jz - 1] &= 0x7fffff; - break; - case 2: - iq[jz - 1] &= 0x3fffff; - break; - } - } - if (ih == 2) { - z = 1.0 - z; - if (carry != 0) - z -= __scalbn(1.0, q0); - } - } - - /* check if recomputation is needed */ - if (z == 0.0) { - j = 0; - for (i = jz - 1; i >= jk; i--) j |= iq[i]; - if (j == 0) { /* need recomputation */ - for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */ - - for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */ - f[jx + i] = (double)ipio2[jv + i]; - for (j = 0, fw = 0.0; j <= jx; j++) - fw += x[j] * f[jx + i - j]; - q[i] = fw; - } - jz += k; - goto recompute; - } - } - - /* chop off zero terms */ - if (z == 0.0) { - jz -= 1; - q0 -= 24; - while (iq[jz] == 0) { - jz--; - q0 -= 24; - } - } else { /* break z into 24-bit if necessary */ - z = __scalbn(z, -q0); - if (z >= 0x1p24) { - fw = (double)(INT32)(0x1p-24 * z); - iq[jz] = (INT32)(z - 0x1p24 * fw); - jz += 1; - q0 += 24; - iq[jz] = (INT32)fw; - } else - iq[jz] = (INT32)z; - } - - /* convert integer "bit" chunk to floating-point value */ - fw = __scalbn(1.0, q0); - for (i = jz; i >= 0; i--) { - q[i] = fw * (double)iq[i]; - fw *= 0x1p-24; - } - - /* compute PIo2[0,...,jp]*q[jz,...,0] */ - for(i = jz; i >= 0; i--) { - for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++) - fw += PIo2[k] * q[i + k]; - fq[jz - i] = fw; - } - - /* compress fq[] into y[] */ - switch(prec) { - case 0: - fw = 0.0; - for (i = jz; i >= 0; i--) - fw += fq[i]; - y[0] = ih == 0 ? fw : -fw; - break; - case 1: - case 2: - fw = 0.0; - for (i = jz; i >= 0; i--) - fw += fq[i]; - fw = (double)fw; - y[0] = ih==0 ? fw : -fw; - fw = fq[0] - fw; - for (i = 1; i <= jz; i++) - fw += fq[i]; - y[1] = ih == 0 ? fw : -fw; - break; - case 3: /* painful */ - for (i = jz; i > 0; i--) { - fw = fq[i - 1] + fq[i]; - fq[i] += fq[i - 1] - fw; - fq[i - 1] = fw; - } - for (i = jz; i > 1; i--) { - fw = fq[i - 1] + fq[i]; - fq[i] += fq[i - 1] - fw; - fq[i - 1] = fw; - } - for (fw = 0.0, i = jz; i >= 2; i--) - fw += fq[i]; - if (ih == 0) { - y[0] = fq[0]; - y[1] = fq[1]; - y[2] = fw; - } else { - y[0] = -fq[0]; - y[1] = -fq[1]; - y[2] = -fw; - } - } - return n & 7; -} - /* Copied from musl: src/math/__rem_pio2.c */ static int __rem_pio2(double x, double *y) { diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c index 239d6005a04..1cec98ae3d4 100644 --- a/dlls/msvcrt/unixlib.c +++ b/dlls/msvcrt/unixlib.c @@ -121,14 +121,6 @@ static float CDECL unix_atanhf(float x) #endif } -/********************************************************************* - * cosf - */ -static float CDECL unix_cosf( float x ) -{ - return cosf( x ); -} - /********************************************************************* * cosh */ @@ -477,7 +469,6 @@ static const struct unix_funcs funcs = unix_asinhf, unix_atanh, unix_atanhf, - unix_cosf, unix_cosh, unix_coshf, unix_exp, diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h index 7cce30885a9..d34117eeb32 100644 --- a/dlls/msvcrt/unixlib.h +++ b/dlls/msvcrt/unixlib.h @@ -29,7 +29,6 @@ struct unix_funcs float (CDECL *asinhf)(float x); double (CDECL *atanh)(double x); float (CDECL *atanhf)(float x); - float (CDECL *cosf)(float x); double (CDECL *cosh)(double x); float (CDECL *coshf)(float x); double (CDECL *exp)(double x);