msvcrt: Import tgamma implementation from musl.

configure

@@ -19623,7 +19623,6 @@ for ac_func in \
 	exp2 \
 	exp2f \
 	fmaf \
-	tgamma \
 	tgammaf
 
 do :

configure.ac

@@ -2662,7 +2662,6 @@ AC_CHECK_FUNCS(\
 	exp2 \
 	exp2f \
 	fmaf \
-	tgamma \
 	tgammaf
 )
 LIBS="$ac_save_LIBS"

dlls/msvcrt/math.c

@@ -8647,16 +8647,15 @@ end:
 /* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */
 static double sin_pi(double x)
 {
-    static const double pi = 3.14159265358979311600e+00;
     int n;
 
     /* spurious inexact if odd int */
     x = 2.0 * (x * 0.5 - floor(x * 0.5)); /* x mod 2.0 */
 
-    n = (int)(x * 4.0);
+    n = x * 4.0;
     n = (n + 1) / 2;
     x -= n * 0.5f;
-    x *= pi;
+    x *= M_PI;
 
     switch (n) {
     default: /* case 4: */
@@ -9047,12 +9046,127 @@ float CDECL lgammaf(float x)
     return r;
 }
 
+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)
 {
-    return unix_funcs->tgamma( 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;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

@@ -111,19 +111,6 @@ static float CDECL unix_powf( float x, float y )
 }
 
 /*********************************************************************
- *      tgamma
- */
-static double CDECL unix_tgamma(double x)
-{
-#ifdef HAVE_TGAMMA
-    return tgamma(x);
-#else
-    FIXME( "not implemented\n" );
-    return 0;
-#endif
-}
-
-/*********************************************************************
  *      tgammaf
  */
 static float CDECL unix_tgammaf(float x)
@@ -145,7 +132,6 @@ static const struct unix_funcs funcs =
     unix_fmaf,
     unix_pow,
     unix_powf,
-    unix_tgamma,
     unix_tgammaf,
 };
 

dlls/msvcrt/unixlib.h

@@ -30,7 +30,6 @@ struct unix_funcs
     float           (CDECL *fmaf)(float x, float y, float z);
     double          (CDECL *pow)(double x, double y);
     float           (CDECL *powf)(float x, float y);
-    double          (CDECL *tgamma)(double x);
     float           (CDECL *tgammaf)(float x);
 };
 

include/config.h.in

@@ -929,9 +929,6 @@
 /* Define to 1 if you have the <termios.h> header file. */
 #undef HAVE_TERMIOS_H
 
-/* Define to 1 if you have the `tgamma' function. */
-#undef HAVE_TGAMMA
-
 /* Define to 1 if you have the `tgammaf' function. */
 #undef HAVE_TGAMMAF