msvcrt: Import expm1 implementation from musl.

Signed-off-by: Piotr Caban <piotr@codeweavers.com>
Signed-off-by: Alexandre Julliard <julliard@winehq.org>

configure

@@ -19618,7 +19618,6 @@ fi
 for ac_func in \
 	exp2 \
 	exp2f \
-	expm1 \
 	expm1f \
 	fma \
 	fmaf \

configure.ac

@@ -2658,7 +2658,6 @@ fi
 AC_CHECK_FUNCS(\
 	exp2 \
 	exp2f \
-	expm1 \
 	expm1f \
 	fma \
 	fmaf \

dlls/msvcrt/math.c

@@ -5758,12 +5758,100 @@ float CDECL exp2f(float x)
 
 /*********************************************************************
  *      expm1 (MSVCR120.@)
+ *
+ * Copied from musl: src/math/expm1.c
  */
 double CDECL expm1(double x)
 {
-    double ret = unix_funcs->expm1( x );
-    if (isfinite(x) && !isfinite(ret)) *_errno() = ERANGE;
-    return ret;
+    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;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

@@ -98,18 +98,6 @@ static float CDECL unix_exp2f( float x )
 #endif
 }
 
-/*********************************************************************
- *      expm1
- */
-static double CDECL unix_expm1(double x)
-{
-#ifdef HAVE_EXPM1
-    return expm1(x);
-#else
-    return exp(x) - 1;
-#endif
-}
-
 /*********************************************************************
  *      expm1f
  */
@@ -374,7 +362,6 @@ static const struct unix_funcs funcs =
     unix_expf,
     unix_exp2,
     unix_exp2f,
-    unix_expm1,
     unix_expm1f,
     unix_fma,
     unix_fmaf,

dlls/msvcrt/unixlib.h

@@ -29,7 +29,6 @@ struct unix_funcs
     float           (CDECL *expf)(float x);
     double          (CDECL *exp2)(double x);
     float           (CDECL *exp2f)(float x);
-    double          (CDECL *expm1)(double x);
     float           (CDECL *expm1f)(float x);
     double          (CDECL *fma)(double x, double y, double z);
     float           (CDECL *fmaf)(float x, float y, float z);

include/config.h.in

@@ -98,9 +98,6 @@
 /* Define to 1 if you have the `exp2f' function. */
 #undef HAVE_EXP2F
 
-/* Define to 1 if you have the `expm1' function. */
-#undef HAVE_EXPM1
-
 /* Define to 1 if you have the `expm1f' function. */
 #undef HAVE_EXPM1F