Module: wine Branch: master Commit: 16fa6b78460c83c08d28aa0a77d6077ea5010650 URL: https://source.winehq.org/git/wine.git/?a=commit;h=16fa6b78460c83c08d28aa0a7...
Author: Piotr Caban piotr@codeweavers.com Date: Tue Jun 8 21:19:27 2021 +0200
msvcrt: Import tgamma implementation from musl.
Signed-off-by: Piotr Caban piotr@codeweavers.com Signed-off-by: Alexandre Julliard julliard@winehq.org
---
configure | 1 - configure.ac | 1 - dlls/msvcrt/math.c | 122 ++++++++++++++++++++++++++++++++++++++++++++++++-- dlls/msvcrt/unixlib.c | 14 ------ dlls/msvcrt/unixlib.h | 1 - include/config.h.in | 3 -- 6 files changed, 118 insertions(+), 24 deletions(-)
diff --git a/configure b/configure index 9acc579ca11..c0d24695456 100755 --- a/configure +++ b/configure @@ -19623,7 +19623,6 @@ for ac_func in \ exp2 \ exp2f \ fmaf \ - tgamma \ tgammaf
do : diff --git a/configure.ac b/configure.ac index bf7efa0690c..748913f7216 100644 --- a/configure.ac +++ b/configure.ac @@ -2662,7 +2662,6 @@ AC_CHECK_FUNCS(\ exp2 \ exp2f \ fmaf \ - tgamma \ tgammaf ) LIBS="$ac_save_LIBS" diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index a471f07f046..d419c34a947 100644 --- a/dlls/msvcrt/math.c +++ b/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; }
/********************************************************************* diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c index 012ec6f6b85..26af9812f88 100644 --- a/dlls/msvcrt/unixlib.c +++ b/dlls/msvcrt/unixlib.c @@ -110,19 +110,6 @@ static float CDECL unix_powf( float x, float y ) return powf( x, y ); }
-/********************************************************************* - * tgamma - */ -static double CDECL unix_tgamma(double x) -{ -#ifdef HAVE_TGAMMA - return tgamma(x); -#else - FIXME( "not implemented\n" ); - return 0; -#endif -} - /********************************************************************* * tgammaf */ @@ -145,7 +132,6 @@ static const struct unix_funcs funcs = unix_fmaf, unix_pow, unix_powf, - unix_tgamma, unix_tgammaf, };
diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h index a9c569d6f63..f119a0a35d0 100644 --- a/dlls/msvcrt/unixlib.h +++ b/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); };
diff --git a/include/config.h.in b/include/config.h.in index 537dda5f68d..aa74e710e06 100644 --- a/include/config.h.in +++ b/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
-
Alexandre Julliard