Module: wine Branch: master Commit: 9a459dc5afa417351de84f9ede4e14aab2c5a590 URL: https://source.winehq.org/git/wine.git/?a=commit;h=9a459dc5afa417351de84f9ed...
Author: Piotr Caban piotr@codeweavers.com Date: Thu Jun 10 19:04:36 2021 +0200
msvcrt: Import expf implementation from musl.
Signed-off-by: Piotr Caban piotr@codeweavers.com Signed-off-by: Alexandre Julliard julliard@winehq.org
---
dlls/msvcrt/math.c | 112 +++++++++++++++++++++++++++++++++++--------------- dlls/msvcrt/unixlib.c | 9 ---- dlls/msvcrt/unixlib.h | 1 - 3 files changed, 78 insertions(+), 44 deletions(-)
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c index ae943e845e8..eaaeb9bec12 100644 --- a/dlls/msvcrt/math.c +++ b/dlls/msvcrt/math.c @@ -668,6 +668,38 @@ static float __cosdf(double x) r = C2 + z * C3; return ((1.0 + z * C0) + w * C1) + (w * z) * r; } + +/* Based on musl implementation: src/math/round.c */ +static double __round(double x) +{ + ULONGLONG llx = *(ULONGLONG*)&x, tmp; + int e = (llx >> 52 & 0x7ff) - 0x3ff; + + if (e >= 52) + return x; + if (e < -1) + return 0 * x; + else if (e == -1) + return signbit(x) ? -1 : 1; + + tmp = 0x000fffffffffffffULL >> e; + if (!(llx & tmp)) + return x; + llx += 0x0008000000000000ULL >> e; + llx &= ~tmp; + return *(double*)&llx; +} + +static const UINT64 exp2f_T[] = { + 0x3ff0000000000000ULL, 0x3fefd9b0d3158574ULL, 0x3fefb5586cf9890fULL, 0x3fef9301d0125b51ULL, + 0x3fef72b83c7d517bULL, 0x3fef54873168b9aaULL, 0x3fef387a6e756238ULL, 0x3fef1e9df51fdee1ULL, + 0x3fef06fe0a31b715ULL, 0x3feef1a7373aa9cbULL, 0x3feedea64c123422ULL, 0x3feece086061892dULL, + 0x3feebfdad5362a27ULL, 0x3feeb42b569d4f82ULL, 0x3feeab07dd485429ULL, 0x3feea47eb03a5585ULL, + 0x3feea09e667f3bcdULL, 0x3fee9f75e8ec5f74ULL, 0x3feea11473eb0187ULL, 0x3feea589994cce13ULL, + 0x3feeace5422aa0dbULL, 0x3feeb737b0cdc5e5ULL, 0x3feec49182a3f090ULL, 0x3feed503b23e255dULL, + 0x3feee89f995ad3adULL, 0x3feeff76f2fb5e47ULL, 0x3fef199bdd85529cULL, 0x3fef3720dcef9069ULL, + 0x3fef5818dcfba487ULL, 0x3fef7c97337b9b5fULL, 0x3fefa4afa2a490daULL, 0x3fefd0765b6e4540ULL +}; #endif
#ifndef __i386__ @@ -1134,11 +1166,50 @@ float CDECL coshf( float x ) */ float CDECL expf( float x ) { - float ret = unix_funcs->expf( x ); - if (isnan(x)) return math_error(_DOMAIN, "expf", x, 0, ret); - if (isfinite(x) && !ret) return math_error(_UNDERFLOW, "expf", x, 0, ret); - if (isfinite(x) && !isfinite(ret)) return math_error(_OVERFLOW, "expf", x, 0, ret); - return ret; + static const double C[] = { + 0x1.c6af84b912394p-5 / (1 << 5) / (1 << 5) / (1 << 5), + 0x1.ebfce50fac4f3p-3 / (1 << 5) / (1 << 5), + 0x1.62e42ff0c52d6p-1 / (1 << 5) + }; + static const double invln2n = 0x1.71547652b82fep+0 * (1 << 5); + + double kd, z, r, r2, y, s; + UINT32 abstop; + UINT64 ki, t; + + abstop = (*(UINT32*)&x >> 20) & 0x7ff; + if (abstop >= 0x42b) { + /* |x| >= 88 or x is nan. */ + if (*(UINT32*)&x == 0xff800000) + return 0.0f; + if (abstop >= 0x7f8) + return x + x; + if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */ + return math_error(_OVERFLOW, "expf", x, 0, x * FLT_MAX); + if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */ + return math_error(_UNDERFLOW, "expf", x, 0, fp_barrierf(FLT_MIN) * FLT_MIN); + } + + /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k. */ + z = invln2n * x; + + /* Round and convert z to int, the result is in [-150*N, 128*N] and + ideally ties-to-even rule is used, otherwise the magnitude of r + can be bigger which gives larger approximation error. */ + kd = __round(z); + ki = kd; + r = z - kd; + + /* exp(x) = 2^(k/N) * 2^(r/N) ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ + t = exp2f_T[ki % (1 << 5)]; + t += ki << (52 - 5); + s = *(double*)&t; + z = C[0] * r + C[1]; + r2 = r * r; + y = C[2] * r + 1; + y = z * r2 + y; + y = y * s; + return y; }
/********************************************************************* @@ -7030,16 +7101,6 @@ double CDECL exp2(double x) */ float CDECL exp2f(float x) { - static const UINT64 T[] = { - 0x3ff0000000000000ULL, 0x3fefd9b0d3158574ULL, 0x3fefb5586cf9890fULL, 0x3fef9301d0125b51ULL, - 0x3fef72b83c7d517bULL, 0x3fef54873168b9aaULL, 0x3fef387a6e756238ULL, 0x3fef1e9df51fdee1ULL, - 0x3fef06fe0a31b715ULL, 0x3feef1a7373aa9cbULL, 0x3feedea64c123422ULL, 0x3feece086061892dULL, - 0x3feebfdad5362a27ULL, 0x3feeb42b569d4f82ULL, 0x3feeab07dd485429ULL, 0x3feea47eb03a5585ULL, - 0x3feea09e667f3bcdULL, 0x3fee9f75e8ec5f74ULL, 0x3feea11473eb0187ULL, 0x3feea589994cce13ULL, - 0x3feeace5422aa0dbULL, 0x3feeb737b0cdc5e5ULL, 0x3feec49182a3f090ULL, 0x3feed503b23e255dULL, - 0x3feee89f995ad3adULL, 0x3feeff76f2fb5e47ULL, 0x3fef199bdd85529cULL, 0x3fef3720dcef9069ULL, - 0x3fef5818dcfba487ULL, 0x3fef7c97337b9b5fULL, 0x3fefa4afa2a490daULL, 0x3fefd0765b6e4540ULL - }; static const double C[] = { 0x1.c6af84b912394p-5, 0x1.ebfce50fac4f3p-3, 0x1.62e42ff0c52d6p-1 }; @@ -7074,7 +7135,7 @@ float CDECL exp2f(float x) r = xd - kd;
/* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */ - t = T[ki % (1 << 5)]; + t = exp2f_T[ki % (1 << 5)]; t += ki << (52 - 5); s = *(double*)&t; z = C[0] * r + C[1]; @@ -7682,27 +7743,10 @@ __int64 CDECL llrintf(float x)
/********************************************************************* * round (MSVCR120.@) - * - * Based on musl implementation: src/math/round.c */ double CDECL round(double x) { - ULONGLONG llx = *(ULONGLONG*)&x, tmp; - int e = (llx >> 52 & 0x7ff) - 0x3ff; - - if (e >= 52) - return x; - if (e < -1) - return 0 * x; - else if (e == -1) - return signbit(x) ? -1 : 1; - - tmp = 0x000fffffffffffffULL >> e; - if (!(llx & tmp)) - return x; - llx += 0x0008000000000000ULL >> e; - llx &= ~tmp; - return *(double*)&llx; + return __round(x); }
/********************************************************************* diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c index 0adf402068f..0498a736812 100644 --- a/dlls/msvcrt/unixlib.c +++ b/dlls/msvcrt/unixlib.c @@ -50,14 +50,6 @@ static double CDECL unix_exp( double x ) return exp( x ); }
-/********************************************************************* - * expf - */ -static float CDECL unix_expf( float x ) -{ - return expf( x ); -} - /********************************************************************* * exp2 */ @@ -89,7 +81,6 @@ static float CDECL unix_powf( float x, float y ) static const struct unix_funcs funcs = { unix_exp, - unix_expf, unix_exp2, unix_pow, unix_powf, diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h index c10f25e73f9..21bb5e6c130 100644 --- a/dlls/msvcrt/unixlib.h +++ b/dlls/msvcrt/unixlib.h @@ -24,7 +24,6 @@ struct unix_funcs { double (CDECL *exp)(double x); - float (CDECL *expf)(float x); double (CDECL *exp2)(double x); double (CDECL *pow)(double x, double y); float (CDECL *powf)(float x, float y);