msvcrt: Use the tgamma()/tgammaf() implementation from the bundled musl library.

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:

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))