Piotr Caban : msvcrt: Import cosh implementation from musl.
Alexandre Julliard
julliard at winehq.org
Mon May 31 15:58:35 CDT 2021
Module: wine
Branch: master
Commit: b431d7515da309ff85e3677fe98fce2218c890ab
URL: https://source.winehq.org/git/wine.git/?a=commit;h=b431d7515da309ff85e3677fe98fce2218c890ab
Author: Piotr Caban <piotr at codeweavers.com>
Date: Sat May 29 19:13:39 2021 +0200
msvcrt: Import cosh implementation from musl.
Signed-off-by: Piotr Caban <piotr at codeweavers.com>
Signed-off-by: Alexandre Julliard <julliard at 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);
More information about the wine-cvs
mailing list