From: Alexandre Julliard julliard@winehq.org
--- dlls/msvcrt/math.c | 155 ------------------------------------ libs/musl/src/math/tgamma.c | 13 ++- 2 files changed, 10 insertions(+), 158 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 20ce8bf9599..0c0e9df41c7 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -8996,161 +8996,6 @@ float CDECL remainderf(float x, float y) return remquof(x, y, &q); }
-/* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ -static double sin_pi(double x) -{ - int n; - - /* spurious inexact if odd int */ - x = 2.0 * (x * 0.5 - floor(x * 0.5)); /* x mod 2.0 */ - - n = x * 4.0; - n = (n + 1) / 2; - x -= n * 0.5f; - x *= M_PI; - - switch (n) { - default: /* case 4: */ - case 0: return __sin(x, 0.0, 0); - case 1: return __cos(x, 0.0); - case 2: return __sin(-x, 0.0, 0); - case 3: return -__cos(x, 0.0); - } -} - -static double tgamma_S(double x) -{ - static const double Snum[] = { - 23531376880.410759688572007674451636754734846804940, - 42919803642.649098768957899047001988850926355848959, - 35711959237.355668049440185451547166705960488635843, - 17921034426.037209699919755754458931112671403265390, - 6039542586.3520280050642916443072979210699388420708, - 1439720407.3117216736632230727949123939715485786772, - 248874557.86205415651146038641322942321632125127801, - 31426415.585400194380614231628318205362874684987640, - 2876370.6289353724412254090516208496135991145378768, - 186056.26539522349504029498971604569928220784236328, - 8071.6720023658162106380029022722506138218516325024, - 210.82427775157934587250973392071336271166969580291, - 2.5066282746310002701649081771338373386264310793408, - }; - static const double Sden[] = { - 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, - 2637558, 357423, 32670, 1925, 66, 1, - }; - - double num = 0, den = 0; - int i; - - /* to avoid overflow handle large x differently */ - if (x < 8) - for (i = ARRAY_SIZE(Snum) - 1; i >= 0; i--) { - num = num * x + Snum[i]; - den = den * x + Sden[i]; - } - else - for (i = 0; i < ARRAY_SIZE(Snum); i++) { - num = num / x + Snum[i]; - den = den / x + Sden[i]; - } - return num / den; -} - -/********************************************************************* - * tgamma (MSVCR120.@) - * - * Copied from musl: src/math/tgamma.c - */ -double CDECL tgamma(double x) -{ - static const double gmhalf = 5.524680040776729583740234375; - static const double fact[] = { - 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0, - 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0, - 355687428096000.0, 6402373705728000.0, 121645100408832000.0, - 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0, - }; - - union {double f; UINT64 i;} u = {x}; - double absx, y, dy, z, r; - UINT32 ix = u.i >> 32 & 0x7fffffff; - int sign = u.i >> 63; - - /* special cases */ - if (ix >= 0x7ff00000) { - /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ - if (u.i == 0xfff0000000000000ULL) - *_errno() = EDOM; - return x + INFINITY; - } - if (ix < (0x3ff - 54) << 20) { - /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ - if (x == 0.0) - *_errno() = ERANGE; - return 1 / x; - } - - /* integer arguments */ - /* raise inexact when non-integer */ - if (x == floor(x)) { - if (sign) { - *_errno() = EDOM; - return 0 / (x - x); - } - if (x <= ARRAY_SIZE(fact)) - return fact[(int)x - 1]; - } - - /* x >= 172: tgamma(x)=inf with overflow */ - /* x =< -184: tgamma(x)=+-0 with underflow */ - if (ix >= 0x40670000) { /* |x| >= 184 */ - *_errno() = ERANGE; - if (sign) { - fp_barrierf(0x1p-126 / x); - return 0; - } - x *= 0x1p1023; - return x; - } - - absx = sign ? -x : x; - - /* handle the error of x + g - 0.5 */ - y = absx + gmhalf; - if (absx > gmhalf) { - dy = y - absx; - dy -= gmhalf; - } else { - dy = y - gmhalf; - dy -= absx; - } - - z = absx - 0.5; - r = tgamma_S(absx) * exp(-y); - if (x < 0) { - /* reflection formula for negative x */ - /* sinpi(absx) is not 0, integers are already handled */ - r = -M_PI / (sin_pi(absx) * absx * r); - dy = -dy; - z = -z; - } - r += dy * (gmhalf + 0.5) * r / y; - z = pow(y, 0.5 * z); - y = r * z * z; - return y; -} - -/********************************************************************* - * tgammaf (MSVCR120.@) - * - * Copied from musl: src/math/tgammaf.c - */ -float CDECL tgammaf(float x) -{ - return tgamma(x); -} - /********************************************************************* * _except1 (MSVCR120.@) * TODO: diff --git a/libs/musl/src/math/tgamma.c b/libs/musl/src/math/tgamma.c index 94c8d14f000..8bea688565b 100644 --- a/libs/musl/src/math/tgamma.c +++ b/libs/musl/src/math/tgamma.c @@ -114,18 +114,24 @@ double __cdecl tgamma(double x) int sign = u.i>>63;
/* special cases */ - if (ix >= 0x7ff00000) + if (ix >= 0x7ff00000) { /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ + if (u.i == 0xfff0000000000000ULL) errno = EDOM; return x + INFINITY; - if (ix < (0x3ff-54)<<20) + } + if (ix < (0x3ff-54)<<20) { /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ + if (x == 0.0) errno = ERANGE; return 1/x; + }
/* integer arguments */ /* raise inexact when non-integer */ if (x == floor(x)) { - if (sign) + if (sign) { + errno = EDOM; return 0/0.0; + } if (x <= sizeof fact/sizeof *fact) return fact[(int)x - 1]; } @@ -133,6 +139,7 @@ double __cdecl tgamma(double x) /* x >= 172: tgamma(x)=inf with overflow */ /* x =< -184: tgamma(x)=+-0 with underflow */ if (ix >= 0x40670000) { /* |x| >= 184 */ + errno = ERANGE; if (sign) { FORCE_EVAL((float)(0x1p-126/x)); if (floor(x) * 0.5 == floor(x * 0.5))