Piotr Caban : msvcrt: Import log1p implementation from musl.

[Alexandre Julliard]

Module: wine
Branch: master
Commit: 2e6fa0a49804da29e683c72826960a95651ef41e
URL:    https://source.winehq.org/git/wine.git/?a=commit;h=2e6fa0a49804da29e683c72826960a95651ef41e

Author: Piotr Caban <piotr at codeweavers.com>
Date:   Fri Jun  4 18:25:26 2021 +0200

msvcrt: Import log1p implementation from musl.

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

---

 configure             |  1 -
 configure.ac          |  1 -
 dlls/msvcrt/math.c    | 70 ++++++++++++++++++++++++++++++++++++++++++++++++---
 dlls/msvcrt/unixlib.c | 13 ----------
 dlls/msvcrt/unixlib.h |  1 -
 include/config.h.in   |  3 ---
 6 files changed, 67 insertions(+), 22 deletions(-)

diff --git a/configure b/configure
index 512ced158fa..e91b319c705 100755
--- a/configure
+++ b/configure
@@ -19622,7 +19622,6 @@ for ac_func in \
 	fmaf \
 	lgamma \
 	lgammaf \
-	log1p \
 	log1pf \
 	log2 \
 	log2f \
diff --git a/configure.ac b/configure.ac
index 23075b36a42..accdbaf6216 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2661,7 +2661,6 @@ AC_CHECK_FUNCS(\
 	fmaf \
 	lgamma \
 	lgammaf \
-	log1p \
 	log1pf \
 	log2 \
 	log2f \
diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c
index 6e6d80f50a0..fbc97aedce4 100644
--- a/dlls/msvcrt/math.c
+++ b/dlls/msvcrt/math.c
@@ -6898,12 +6898,76 @@ float CDECL expm1f(float x)
 
 /*********************************************************************
  *      log1p (MSVCR120.@)
+ *
+ * Copied from musl: src/math/log1p.c
  */
 double CDECL log1p(double x)
 {
-    if (x < -1) *_errno() = EDOM;
-    else if (x == -1) *_errno() = ERANGE;
-    return unix_funcs->log1p( x );
+    static const double ln2_hi = 6.93147180369123816490e-01,
+        ln2_lo = 1.90821492927058770002e-10,
+        Lg1 = 6.666666666666735130e-01,
+        Lg2 = 3.999999999940941908e-01,
+        Lg3 = 2.857142874366239149e-01,
+        Lg4 = 2.222219843214978396e-01,
+        Lg5 = 1.818357216161805012e-01,
+        Lg6 = 1.531383769920937332e-01,
+        Lg7 = 1.479819860511658591e-01;
+
+    union {double f; UINT64 i;} u = {x};
+    double hfsq, f, c, s, z, R, w, t1, t2, dk;
+    UINT32 hx, hu;
+    int k;
+
+    hx = u.i >> 32;
+    k = 1;
+    if (hx < 0x3fda827a || hx >> 31) { /* 1+x < sqrt(2)+ */
+        if (hx >= 0xbff00000) { /* x <= -1.0 */
+            if (x == -1) {
+                *_errno() = ERANGE;
+                return x / 0.0; /* og1p(-1) = -inf */
+            }
+            *_errno() = EDOM;
+            return (x-x) / 0.0; /* log1p(x<-1) = NaN */
+        }
+        if (hx << 1 < 0x3ca00000 << 1) { /* |x| < 2**-53 */
+            fp_barrier(x + 0x1p120f);
+            /* underflow if subnormal */
+            if ((hx & 0x7ff00000) == 0)
+                fp_barrierf(x);
+            return x;
+        }
+        if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+            k = 0;
+            c = 0;
+            f = x;
+        }
+    } else if (hx >= 0x7ff00000)
+        return x;
+    if (k) {
+        u.f = 1 + x;
+        hu = u.i >> 32;
+        hu += 0x3ff00000 - 0x3fe6a09e;
+        k = (int)(hu >> 20) - 0x3ff;
+        /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
+        if (k < 54) {
+            c = k >= 2 ? 1 - (u.f - x) : x - (u.f - 1);
+            c /= u.f;
+        } else
+            c = 0;
+        /* reduce u into [sqrt(2)/2, sqrt(2)] */
+        hu = (hu & 0x000fffff) + 0x3fe6a09e;
+        u.i = (UINT64)hu << 32 | (u.i & 0xffffffff);
+        f = u.f - 1;
+    }
+    hfsq = 0.5 * f * f;
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+    t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+    R = t2 + t1;
+    dk = k;
+    return s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi;
 }
 
 /*********************************************************************
diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c
index c880bcda8e2..f98d08e4909 100644
--- a/dlls/msvcrt/unixlib.c
+++ b/dlls/msvcrt/unixlib.c
@@ -136,18 +136,6 @@ static float CDECL unix_log10f( float x )
     return log10f( x );
 }
 
-/*********************************************************************
- *      log1p
- */
-static double CDECL unix_log1p(double x)
-{
-#ifdef HAVE_LOG1P
-    return log1p(x);
-#else
-    return log(1 + x);
-#endif
-}
-
 /*********************************************************************
  *      log1pf
  */
@@ -237,7 +225,6 @@ static const struct unix_funcs funcs =
     unix_lgammaf,
     unix_log10,
     unix_log10f,
-    unix_log1p,
     unix_log1pf,
     unix_log2,
     unix_log2f,
diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h
index bbec27c1061..285deb22bc2 100644
--- a/dlls/msvcrt/unixlib.h
+++ b/dlls/msvcrt/unixlib.h
@@ -32,7 +32,6 @@ struct unix_funcs
     float           (CDECL *lgammaf)(float x);
     double          (CDECL *log10)(double x);
     float           (CDECL *log10f)(float x);
-    double          (CDECL *log1p)(double x);
     float           (CDECL *log1pf)(float x);
     double          (CDECL *log2)(double x);
     float           (CDECL *log2f)(float x);
diff --git a/include/config.h.in b/include/config.h.in
index 3d1e5e6ca5a..7e6d36112d3 100644
--- a/include/config.h.in
+++ b/include/config.h.in
@@ -408,9 +408,6 @@
 /* Define to 1 if you have the <linux/videodev2.h> header file. */
 #undef HAVE_LINUX_VIDEODEV2_H
 
-/* Define to 1 if you have the `log1p' function. */
-#undef HAVE_LOG1P
-
 /* Define to 1 if you have the `log1pf' function. */
 #undef HAVE_LOG1PF