Module: wine Branch: master Commit: b431d7515da309ff85e3677fe98fce2218c890ab URL: https://source.winehq.org/git/wine.git/?a=commit;h=b431d7515da309ff85e3677fe...
Author: Piotr Caban piotr@codeweavers.com Date: Sat May 29 19:13:39 2021 +0200
msvcrt: Import cosh implementation from musl.
Signed-off-by: Piotr Caban piotr@codeweavers.com Signed-off-by: Alexandre Julliard julliard@winehq.org
---
dlls/msvcrt/math.c | 231 ++++++++++++++++++++++++++++++-------------------- dlls/msvcrt/unixlib.c | 9 -- dlls/msvcrt/unixlib.h | 1 - 3 files changed, 137 insertions(+), 104 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index 06f1dc26838..760f3834761 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -2095,14 +2095,147 @@ double CDECL cos( double x ) } }
+/* Copied from musl: src/math/expm1.c */ +static double CDECL __expm1(double x) +{ + static const double o_threshold = 7.09782712893383973096e+02, + ln2_hi = 6.93147180369123816490e-01, + ln2_lo = 1.90821492927058770002e-10, + invln2 = 1.44269504088896338700e+00, + Q1 = -3.33333333333331316428e-02, + Q2 = 1.58730158725481460165e-03, + Q3 = -7.93650757867487942473e-05, + Q4 = 4.00821782732936239552e-06, + Q5 = -2.01099218183624371326e-07; + + double y, hi, lo, c, t, e, hxs, hfx, r1, twopk; + union {double f; UINT64 i;} u = {x}; + UINT32 hx = u.i >> 32 & 0x7fffffff; + int k, sign = u.i >> 63; + + /* filter out huge and non-finite argument */ + if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ + if (isnan(x)) + return x; + if (isinf(x)) + return sign ? -1 : x; + if (sign) + return math_error(_UNDERFLOW, "exp", x, 0, -1); + if (x > o_threshold) + return math_error(_OVERFLOW, "exp", x, 0, x * 0x1p1023); + } + + /* argument reduction */ + if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ + if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ + if (!sign) { + hi = x - ln2_hi; + lo = ln2_lo; + k = 1; + } else { + hi = x + ln2_hi; + lo = -ln2_lo; + k = -1; + } + } else { + k = invln2 * x + (sign ? -0.5 : 0.5); + t = k; + hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ + lo = t * ln2_lo; + } + x = hi - lo; + c = (hi - x) - lo; + } else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */ + fp_barrier(x + 0x1p120f); + if (hx < 0x00100000) + fp_barrier((float)x); + return x; + } else + k = 0; + + /* x is now in primary range */ + hfx = 0.5 * x; + hxs = x * hfx; + r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); + t = 3.0 - r1 * hfx; + e = hxs * ((r1 - t) / (6.0 - x * t)); + if (k == 0) /* c is 0 */ + return x - (x * e - hxs); + e = x * (e - c) - c; + e -= hxs; + /* exp(x) ~ 2^k (x_reduced - e + 1) */ + if (k == -1) + return 0.5 * (x - e) - 0.5; + if (k == 1) { + if (x < -0.25) + return -2.0 * (e - (x + 0.5)); + return 1.0 + 2.0 * (x - e); + } + u.i = (UINT64)(0x3ff + k) << 52; /* 2^k */ + twopk = u.f; + if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ + y = x - e + 1.0; + if (k == 1024) + y = y * 2.0 * 0x1p1023; + else + y = y * twopk; + return y - 1.0; + } + u.i = (UINT64)(0x3ff - k) << 52; /* 2^-k */ + if (k < 20) + y = (x - e + (1 - u.f)) * twopk; + else + y = (x - (e + u.f) + 1) * twopk; + return y; +} + +static double __expo2(double x, double sign) +{ + static const int k = 2043; + static const double kln2 = 0x1.62066151add8bp+10; + double scale; + + *(UINT64*)&scale = (UINT64)(0x3ff + k / 2) << 52; + return exp(x - kln2) * (sign * scale) * scale; +} + /********************************************************************* * cosh (MSVCRT.@) + * + * Copied from musl: src/math/cosh.c */ double CDECL cosh( double x ) { - double ret = unix_funcs->cosh( x ); - if (isnan(x)) return math_error(_DOMAIN, "cosh", x, 0, ret); - return ret; + UINT64 ux = *(UINT64*)&x; + UINT32 w; + double t; + + /* |x| */ + ux &= (uint64_t)-1 / 2; + x = *(double*)&ux; + w = ux >> 32; + + /* |x| < log(2) */ + if (w < 0x3fe62e42) { + if (w < 0x3ff00000 - (26 << 20)) { + fp_barrier(x + 0x1p120f); + return 1; + } + t = __expm1(x); + return 1 + t * t / (2 * (1 + t)); + } + + /* |x| < log(DBL_MAX) */ + if (w < 0x40862e42) { + t = exp(x); + /* note: if x>log(0x1p26) then the 1/t is not needed */ + return 0.5 * (t + 1 / t); + } + + /* |x| > log(DBL_MAX) or nan */ + /* note: the result is stored to handle overflow */ + t = __expo2(x, 1.0); + return t; }
/********************************************************************* @@ -5758,100 +5891,10 @@ float CDECL exp2f(float x)
/********************************************************************* * expm1 (MSVCR120.@) - * - * Copied from musl: src/math/expm1.c */ double CDECL expm1(double x) { - static const double o_threshold = 7.09782712893383973096e+02, - ln2_hi = 6.93147180369123816490e-01, - ln2_lo = 1.90821492927058770002e-10, - invln2 = 1.44269504088896338700e+00, - Q1 = -3.33333333333331316428e-02, - Q2 = 1.58730158725481460165e-03, - Q3 = -7.93650757867487942473e-05, - Q4 = 4.00821782732936239552e-06, - Q5 = -2.01099218183624371326e-07; - - double y, hi, lo, c, t, e, hxs, hfx, r1, twopk; - union {double f; UINT64 i;} u = {x}; - UINT32 hx = u.i >> 32 & 0x7fffffff; - int k, sign = u.i >> 63; - - /* filter out huge and non-finite argument */ - if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ - if (isnan(x)) - return x; - if (isinf(x)) - return sign ? -1 : x; - if (sign) - return math_error(_UNDERFLOW, "exp", x, 0, -1); - if (x > o_threshold) - return math_error(_OVERFLOW, "exp", x, 0, x * 0x1p1023); - } - - /* argument reduction */ - if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ - if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ - if (!sign) { - hi = x - ln2_hi; - lo = ln2_lo; - k = 1; - } else { - hi = x + ln2_hi; - lo = -ln2_lo; - k = -1; - } - } else { - k = invln2 * x + (sign ? -0.5 : 0.5); - t = k; - hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ - lo = t * ln2_lo; - } - x = hi - lo; - c = (hi - x) - lo; - } else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */ - fp_barrier(x + 0x1p120f); - if (hx < 0x00100000) - fp_barrier((float)x); - return x; - } else - k = 0; - - /* x is now in primary range */ - hfx = 0.5 * x; - hxs = x * hfx; - r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); - t = 3.0 - r1 * hfx; - e = hxs * ((r1 - t) / (6.0 - x * t)); - if (k == 0) /* c is 0 */ - return x - (x * e - hxs); - e = x * (e - c) - c; - e -= hxs; - /* exp(x) ~ 2^k (x_reduced - e + 1) */ - if (k == -1) - return 0.5 * (x - e) - 0.5; - if (k == 1) { - if (x < -0.25) - return -2.0 * (e - (x + 0.5)); - return 1.0 + 2.0 * (x - e); - } - u.i = (UINT64)(0x3ff + k) << 52; /* 2^k */ - twopk = u.f; - if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ - y = x - e + 1.0; - if (k == 1024) - y = y * 2.0 * 0x1p1023; - else - y = y * twopk; - return y - 1.0; - } - u.i = (UINT64)(0x3ff - k) << 52; /* 2^-k */ - if (k < 20) - y = (x - e + (1 - u.f)) * twopk; - else - y = (x - (e + u.f) + 1) * twopk; - return y; + return __expm1(x); }
/********************************************************************* diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c index c0c3a2b1636..6af252e5e9b 100644 --- a/dlls/msvcrt/unixlib.c +++ b/dlls/msvcrt/unixlib.c @@ -42,14 +42,6 @@
WINE_DEFAULT_DEBUG_CHANNEL(msvcrt);
-/********************************************************************* - * cosh - */ -static double CDECL unix_cosh( double x ) -{ - return cosh( x ); -} - /********************************************************************* * coshf */ @@ -344,7 +336,6 @@ static float CDECL unix_tgammaf(float x)
static const struct unix_funcs funcs = { - unix_cosh, unix_coshf, unix_exp, unix_expf, diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h index 6f8347b066a..67dd1f98b84 100644 --- a/dlls/msvcrt/unixlib.h +++ b/dlls/msvcrt/unixlib.h @@ -23,7 +23,6 @@
struct unix_funcs { - double (CDECL *cosh)(double x); float (CDECL *coshf)(float x); double (CDECL *exp)(double x); float (CDECL *expf)(float x);