From b9b19ea2db2266b3cb35d20882e286dc117d04c5 Mon Sep 17 00:00:00 2001 From: Alexandre Julliard Date: Mon, 3 Apr 2023 17:36:09 +0200 Subject: [PATCH] msvcrt: Use the __rem_pio2()/__rem_pio2f() implementation from the bundled musl library. --- dlls/msvcrt/math.c | 424 +------------------------------ libs/musl/src/math/__rem_pio2.c | 2 +- libs/musl/src/math/__rem_pio2f.c | 2 +- 3 files changed, 10 insertions(+), 418 deletions(-) diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 68218b7007b..befd05e9ad0 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -218,212 +218,6 @@ float CDECL _chgsignf( float num ) #endif -/* 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__ /* Copied from musl: src/math/__sindf.c */ static float __sindf(double x) @@ -731,56 +525,7 @@ float CDECL atanf( float x ) return sign ? -z : z; } -/* 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; -} +extern int __rem_pio2f(float x, double *y); /********************************************************************* * cosf (MSVCRT.@) @@ -1627,8 +1372,12 @@ double CDECL atan( double x ) return sign ? -z : z; } -/* Copied from musl: src/math/rint.c */ -static double __rint(double x) +/********************************************************************* + * rint (MSVCR120.@) + * + * Copied from musl: src/math/rint.c + */ +double CDECL rint(double x) { static const double toint = 1 / DBL_EPSILON; @@ -1654,156 +1403,7 @@ static double __rint(double x) return y; } -/* Copied from musl: src/math/__rem_pio2.c */ -static int __rem_pio2(double x, double *y) -{ - static const double pio4 = 0x1.921fb54442d18p-1, - invpio2 = 6.36619772367581382433e-01, - pio2_1 = 1.57079632673412561417e+00, - pio2_1t = 6.07710050650619224932e-11, - pio2_2 = 6.07710050630396597660e-11, - pio2_2t = 2.02226624879595063154e-21, - pio2_3 = 2.02226624871116645580e-21, - pio2_3t = 8.47842766036889956997e-32; - - union {double f; UINT64 i;} u = {x}; - double z, w, t, r, fn, tx[3], ty[2]; - UINT32 ix; - int sign, n, ex, ey, i; - - sign = u.i >> 63; - ix = u.i >> 32 & 0x7fffffff; - if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ - if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ - goto medium; /* cancellation -- use medium case */ - if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ - if (!sign) { - z = x - pio2_1; /* one round good to 85 bits */ - y[0] = z - pio2_1t; - y[1] = (z - y[0]) - pio2_1t; - return 1; - } else { - z = x + pio2_1; - y[0] = z + pio2_1t; - y[1] = (z - y[0]) + pio2_1t; - return -1; - } - } else { - if (!sign) { - z = x - 2 * pio2_1; - y[0] = z - 2 * pio2_1t; - y[1] = (z - y[0]) - 2 * pio2_1t; - return 2; - } else { - z = x + 2 * pio2_1; - y[0] = z + 2 * pio2_1t; - y[1] = (z - y[0]) + 2 * pio2_1t; - return -2; - } - } - } - if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ - if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ - if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ - goto medium; - if (!sign) { - z = x - 3 * pio2_1; - y[0] = z - 3 * pio2_1t; - y[1] = (z - y[0]) - 3 * pio2_1t; - return 3; - } else { - z = x + 3 * pio2_1; - y[0] = z + 3 * pio2_1t; - y[1] = (z - y[0]) + 3 * pio2_1t; - return -3; - } - } else { - if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ - goto medium; - if (!sign) { - z = x - 4 * pio2_1; - y[0] = z - 4 * pio2_1t; - y[1] = (z - y[0]) - 4 * pio2_1t; - return 4; - } else { - z = x + 4 * pio2_1; - y[0] = z + 4 * pio2_1t; - y[1] = (z - y[0]) + 4 * pio2_1t; - return -4; - } - } - } - if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ -medium: - fn = __rint(x * invpio2); - n = (INT32)fn; - r = x - fn * pio2_1; - w = fn * pio2_1t; /* 1st round, good to 85 bits */ - /* Matters with directed rounding. */ - if (r - w < -pio4) { - n--; - fn--; - r = x - fn * pio2_1; - w = fn * pio2_1t; - } else if (r - w > pio4) { - n++; - fn++; - r = x - fn * pio2_1; - w = fn * pio2_1t; - } - y[0] = r - w; - u.f = y[0]; - ey = u.i >> 52 & 0x7ff; - ex = ix >> 20; - if (ex - ey > 16) { /* 2nd round, good to 118 bits */ - t = r; - w = fn * pio2_2; - r = t - w; - w = fn * pio2_2t - ((t - r) - w); - y[0] = r - w; - u.f = y[0]; - ey = u.i >> 52 & 0x7ff; - if (ex - ey > 49) { /* 3rd round, good to 151 bits, covers all cases */ - t = r; - w = fn * pio2_3; - r = t - w; - w = fn * pio2_3t - ((t - r) - w); - y[0] = r - w; - } - } - y[1] = (r - y[0]) - w; - return n; - } - /* - * all other (large) arguments - */ - if (ix >= 0x7ff00000) { /* x is inf or NaN */ - y[0] = y[1] = x - x; - return 0; - } - /* set z = scalbn(|x|,-ilogb(x)+23) */ - u.f = x; - u.i &= (UINT64)-1 >> 12; - u.i |= (UINT64)(0x3ff + 23) << 52; - z = u.f; - for (i = 0; i < 2; i++) { - tx[i] = (double)(INT32)z; - z = (z - tx[i]) * 0x1p24; - } - tx[i] = z; - /* skip zero terms, first term is non-zero */ - while (tx[i] == 0.0) - i--; - n = __rem_pio2_large(tx, ty, (int)(ix >> 20) - (0x3ff + 23), i + 1, 1); - if (sign) { - y[0] = -ty[0]; - y[1] = -ty[1]; - return -n; - } - y[0] = ty[0]; - y[1] = ty[1]; - return n; -} +extern int __rem_pio2(double x, double *y); /* Copied from musl: src/math/__sin.c */ static double __sin(double x, double y, int iy) @@ -4962,14 +4562,6 @@ void __cdecl __libm_sse2_sqrt_precise(void) #if _MSVCR_VER>=120 -/********************************************************************* - * rint (MSVCR120.@) - */ -double CDECL rint(double x) -{ - return __rint(x); -} - /********************************************************************* * rintf (MSVCR120.@) * diff --git a/libs/musl/src/math/__rem_pio2.c b/libs/musl/src/math/__rem_pio2.c index dcf672fbd70..3addef65b35 100644 --- a/libs/musl/src/math/__rem_pio2.c +++ b/libs/musl/src/math/__rem_pio2.c @@ -119,7 +119,7 @@ int __rem_pio2(double x, double *y) if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ medium: /* rint(x/(pi/2)) */ - fn = (double_t)x*invpio2 + toint - toint; + fn = rint(x * invpio2); n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ diff --git a/libs/musl/src/math/__rem_pio2f.c b/libs/musl/src/math/__rem_pio2f.c index e67656431a8..e36fb47f75c 100644 --- a/libs/musl/src/math/__rem_pio2f.c +++ b/libs/musl/src/math/__rem_pio2f.c @@ -52,7 +52,7 @@ int __rem_pio2f(float x, double *y) /* 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 = (double_t)x*invpio2 + toint - toint; + fn = fp_barrier(x * invpio2 + toint) - toint; n = (int32_t)fn; *y = x - fn*pio2_1 - fn*pio2_1t; /* Matters with directed rounding. */