From 7b0e5b2b69cd4a3546d5a95bb77d89de54b83c6e Mon Sep 17 00:00:00 2001 From: Piotr Caban Date: Wed, 19 May 2021 15:26:33 +0200 Subject: [PATCH] msvcrt: Import erfc implementation from musl. Signed-off-by: Piotr Caban Signed-off-by: Alexandre Julliard --- configure | 1 - configure.ac | 1 - dlls/msvcrt/math.c | 128 +++++++++++++++++++++++++++++++++++++++++- dlls/msvcrt/unixlib.c | 13 ----- dlls/msvcrt/unixlib.h | 1 - include/config.h.in | 3 - 6 files changed, 127 insertions(+), 20 deletions(-) diff --git a/configure b/configure index 9eb6cb7fd61..8e15b71d74e 100755 --- a/configure +++ b/configure @@ -19623,7 +19623,6 @@ for ac_func in \ atanh \ atanhf \ erf \ - erfc \ exp2 \ exp2f \ expm1 \ diff --git a/configure.ac b/configure.ac index d21882c7813..90522c22867 100644 --- a/configure.ac +++ b/configure.ac @@ -2663,7 +2663,6 @@ AC_CHECK_FUNCS(\ atanh \ atanhf \ erf \ - erfc \ exp2 \ exp2f \ expm1 \ diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index e0c28e25638..787bb18ddad 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -5215,6 +5215,90 @@ short CDECL _fdtest(float *x) return _fdclass(*x); } +static double erfc1(double x) +{ + static const double erx = 8.45062911510467529297e-01, + pa0 = -2.36211856075265944077e-03, + pa1 = 4.14856118683748331666e-01, + pa2 = -3.72207876035701323847e-01, + pa3 = 3.18346619901161753674e-01, + pa4 = -1.10894694282396677476e-01, + pa5 = 3.54783043256182359371e-02, + pa6 = -2.16637559486879084300e-03, + qa1 = 1.06420880400844228286e-01, + qa2 = 5.40397917702171048937e-01, + qa3 = 7.18286544141962662868e-02, + qa4 = 1.26171219808761642112e-01, + qa5 = 1.36370839120290507362e-02, + qa6 = 1.19844998467991074170e-02; + + double s, P, Q; + + s = fabs(x) - 1; + P = pa0 + s * (pa1 + s * (pa2 + s * (pa3 + s * (pa4 + s * (pa5 + s * pa6))))); + Q = 1 + s * (qa1 + s * (qa2 + s * (qa3 + s * (qa4 + s * (qa5 + s * qa6))))); + return 1 - erx - P / Q; +} + +static double erfc2(uint32_t ix, double x) +{ + static const double ra0 = -9.86494403484714822705e-03, + ra1 = -6.93858572707181764372e-01, + ra2 = -1.05586262253232909814e+01, + ra3 = -6.23753324503260060396e+01, + ra4 = -1.62396669462573470355e+02, + ra5 = -1.84605092906711035994e+02, + ra6 = -8.12874355063065934246e+01, + ra7 = -9.81432934416914548592e+00, + sa1 = 1.96512716674392571292e+01, + sa2 = 1.37657754143519042600e+02, + sa3 = 4.34565877475229228821e+02, + sa4 = 6.45387271733267880336e+02, + sa5 = 4.29008140027567833386e+02, + sa6 = 1.08635005541779435134e+02, + sa7 = 6.57024977031928170135e+00, + sa8 = -6.04244152148580987438e-02, + rb0 = -9.86494292470009928597e-03, + rb1 = -7.99283237680523006574e-01, + rb2 = -1.77579549177547519889e+01, + rb3 = -1.60636384855821916062e+02, + rb4 = -6.37566443368389627722e+02, + rb5 = -1.02509513161107724954e+03, + rb6 = -4.83519191608651397019e+02, + sb1 = 3.03380607434824582924e+01, + sb2 = 3.25792512996573918826e+02, + sb3 = 1.53672958608443695994e+03, + sb4 = 3.19985821950859553908e+03, + sb5 = 2.55305040643316442583e+03, + sb6 = 4.74528541206955367215e+02, + sb7 = -2.24409524465858183362e+01; + + double s, R, S, z; + UINT64 iz; + + if (ix < 0x3ff40000) /* |x| < 1.25 */ + return erfc1(x); + + x = fabs(x); + s = 1 / (x * x); + if (ix < 0x4006db6d) { /* |x| < 1/.35 ~ 2.85714 */ + R = ra0 + s * (ra1 + s * (ra2 + s * (ra3 + s * (ra4 + s * + (ra5 + s * (ra6 + s * ra7)))))); + S = 1.0 + s * (sa1 + s * (sa2 + s * (sa3 + s * (sa4 + s * + (sa5 + s * (sa6 + s * (sa7 + s * sa8))))))); + } else { /* |x| > 1/.35 */ + R = rb0 + s * (rb1 + s * (rb2 + s * (rb3 + s * (rb4 + s * + (rb5 + s * rb6))))); + S = 1.0 + s * (sb1 + s * (sb2 + s * (sb3 + s * (sb4 + s * + (sb5 + s * (sb6 + s * sb7)))))); + } + z = x; + iz = *(ULONGLONG*)&z; + iz &= 0xffffffff00000000ULL; + z = *(double*)&iz; + return exp(-z * z - 0.5625) * exp((z - x) * (z + x) + R / S) / x; +} + /********************************************************************* * erf (MSVCR120.@) */ @@ -5354,10 +5438,52 @@ float CDECL erff(float x) /********************************************************************* * erfc (MSVCR120.@) + * + * Copied from musl: src/math/erf.c */ double CDECL erfc(double x) { - return unix_funcs->erfc( x ); + static const double pp0 = 1.28379167095512558561e-01, + pp1 = -3.25042107247001499370e-01, + pp2 = -2.84817495755985104766e-02, + pp3 = -5.77027029648944159157e-03, + pp4 = -2.37630166566501626084e-05, + qq1 = 3.97917223959155352819e-01, + qq2 = 6.50222499887672944485e-02, + qq3 = 5.08130628187576562776e-03, + qq4 = 1.32494738004321644526e-04, + qq5 = -3.96022827877536812320e-06; + + double r, s, z, y; + UINT32 ix; + int sign; + + ix = *(ULONGLONG*)&x >> 32; + sign = ix >> 31; + ix &= 0x7fffffff; + if (ix >= 0x7ff00000) { + /* erfc(nan)=nan, erfc(+-inf)=0,2 */ + return 2 * sign + 1 / x; + } + if (ix < 0x3feb0000) { /* |x| < 0.84375 */ + if (ix < 0x3c700000) /* |x| < 2**-56 */ + return 1.0 - x; + z = x * x; + r = pp0 + z * (pp1 + z * (pp2 + z * (pp3 + z * pp4))); + s = 1.0 + z * (qq1 + z * (qq2 + z * (qq3 + z * (qq4 + z * qq5)))); + y = r / s; + if (sign || ix < 0x3fd00000) { /* x < 1/4 */ + return 1.0 - (x + x * y); + } + return 0.5 - (x - 0.5 + x * y); + } + if (ix < 0x403c0000) { /* 0.84375 <= |x| < 28 */ + return sign ? 2 - erfc2(ix, x) : erfc2(ix, x); + } + if (sign) + return 2 - DBL_MIN; + *_errno() = ERANGE; + return fp_barrier(DBL_MIN) * DBL_MIN; } /********************************************************************* diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c index 9607ca94e7c..59f8ce4cee7 100644 --- a/dlls/msvcrt/unixlib.c +++ b/dlls/msvcrt/unixlib.c @@ -173,18 +173,6 @@ static double CDECL unix_erf(double x) #endif } -/********************************************************************* - * erfc - */ -static double CDECL unix_erfc(double x) -{ -#ifdef HAVE_ERFC - return erfc(x); -#else - return 1 - unix_erf(x); -#endif -} - /********************************************************************* * exp */ @@ -554,7 +542,6 @@ static const struct unix_funcs funcs = unix_cosh, unix_coshf, unix_erf, - unix_erfc, unix_exp, unix_expf, unix_exp2, diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h index 58ec53d285c..9793fca93bd 100644 --- a/dlls/msvcrt/unixlib.h +++ b/dlls/msvcrt/unixlib.h @@ -34,7 +34,6 @@ struct unix_funcs double (CDECL *cosh)(double x); float (CDECL *coshf)(float x); double (CDECL *erf)(double x); - double (CDECL *erfc)(double x); double (CDECL *exp)(double x); float (CDECL *expf)(float x); double (CDECL *exp2)(double x); diff --git a/include/config.h.in b/include/config.h.in index b0e1fb55332..befb0af9f7e 100644 --- a/include/config.h.in +++ b/include/config.h.in @@ -113,9 +113,6 @@ /* Define to 1 if you have the `erf' function. */ #undef HAVE_ERF -/* Define to 1 if you have the `erfc' function. */ -#undef HAVE_ERFC - /* Define to 1 if you have the `exp2' function. */ #undef HAVE_EXP2