Piotr Caban : msvcrt: Import expf implementation from musl.
Alexandre Julliard
julliard at winehq.org
Thu Jun 10 16:04:51 CDT 2021
Module: wine
Branch: master
Commit: 9a459dc5afa417351de84f9ede4e14aab2c5a590
URL: https://source.winehq.org/git/wine.git/?a=commit;h=9a459dc5afa417351de84f9ede4e14aab2c5a590
Author: Piotr Caban <piotr at codeweavers.com>
Date: Thu Jun 10 19:04:36 2021 +0200
msvcrt: Import expf 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 | 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);
More information about the wine-cvs
mailing list