Skip to content

W
wine-cw
Commit b431d751 authored by Piotr Caban's avatar Piotr Caban Committed by Alexandre Julliard
Browse files
Options

msvcrt: Import cosh implementation from musl.

Signed-off-by: 's avatarPiotr Caban <piotr@codeweavers.com> Signed-off-by: 's avatarAlexandre Julliard <julliard@winehq.org>
parent 7d0d7690
Hide whitespace changes
Inline Side-by-side
Showing with 137 additions and 104 deletions
dlls/msvcrt/math.c
View file @ b431d751
......@@ -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);
}
/*********************************************************************
......
dlls/msvcrt/unixlib.c
View file @ b431d751
......@@ -43,14 +43,6 @@
WINE_DEFAULT_DEBUG_CHANNEL(msvcrt);
/*********************************************************************
* cosh
*/
static double CDECL unix_cosh( double x )
{
return cosh( x );
}
/*********************************************************************
* coshf
*/
static float CDECL unix_coshf( float x )
......@@ -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,
......
dlls/msvcrt/unixlib.h
View file @ b431d751
......@@ -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);
......
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