Commit 7b0e5b2b authored by Piotr Caban's avatar Piotr Caban Committed by Alexandre Julliard

msvcrt: Import erfc implementation from musl.

parent 98525b04
......@@ -19623,7 +19623,6 @@ for ac_func in \
atanh \
atanhf \
erf \
erfc \
exp2 \
exp2f \
expm1 \
......
......@@ -2663,7 +2663,6 @@ AC_CHECK_FUNCS(\
atanh \
atanhf \
erf \
erfc \
exp2 \
exp2f \
expm1 \
......
......@@ -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;
}
/*********************************************************************
......
......@@ -174,18 +174,6 @@ static double CDECL unix_erf(double x)
}
/*********************************************************************
* erfc
*/
static double CDECL unix_erfc(double x)
{
#ifdef HAVE_ERFC
return erfc(x);
#else
return 1 - unix_erf(x);
#endif
}
/*********************************************************************
* exp
*/
static double CDECL unix_exp( double x )
......@@ -554,7 +542,6 @@ static const struct unix_funcs funcs =
unix_cosh,
unix_coshf,
unix_erf,
unix_erfc,
unix_exp,
unix_expf,
unix_exp2,
......
......@@ -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);
......
......@@ -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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment