msvcrt: Import logf implementation from musl.

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

dlls/msvcrt/math.c

@@ -1211,13 +1211,83 @@ float CDECL fmodf( float x, float y )
 
 /*********************************************************************
  *      logf (MSVCRT.@)
+ *
+ * Copied from musl: src/math/logf.c src/math/logf_data.c
  */
 float CDECL logf( float x )
 {
-    float ret = unix_funcs->logf( x );
-    if (x < 0.0) return math_error(_DOMAIN, "logf", x, 0, ret);
-    if (x == 0.0) return math_error(_SING, "logf", x, 0, ret);
-  return ret;
+    static const double Ln2 = 0x1.62e42fefa39efp-1;
+    static const double A[] = {
+        -0x1.00ea348b88334p-2,
+        0x1.5575b0be00b6ap-2,
+        -0x1.ffffef20a4123p-2
+    };
+    static const struct {
+        double invc, logc;
+    } T[] = {
+        { 0x1.661ec79f8f3bep+0, -0x1.57bf7808caadep-2 },
+        { 0x1.571ed4aaf883dp+0, -0x1.2bef0a7c06ddbp-2 },
+        { 0x1.49539f0f010bp+0, -0x1.01eae7f513a67p-2 },
+        { 0x1.3c995b0b80385p+0, -0x1.b31d8a68224e9p-3 },
+        { 0x1.30d190c8864a5p+0, -0x1.6574f0ac07758p-3 },
+        { 0x1.25e227b0b8eap+0, -0x1.1aa2bc79c81p-3 },
+        { 0x1.1bb4a4a1a343fp+0, -0x1.a4e76ce8c0e5ep-4 },
+        { 0x1.12358f08ae5bap+0, -0x1.1973c5a611cccp-4 },
+        { 0x1.0953f419900a7p+0, -0x1.252f438e10c1ep-5 },
+        { 0x1p+0, 0x0p+0 },
+        { 0x1.e608cfd9a47acp-1, 0x1.aa5aa5df25984p-5 },
+        { 0x1.ca4b31f026aap-1, 0x1.c5e53aa362eb4p-4 },
+        { 0x1.b2036576afce6p-1, 0x1.526e57720db08p-3 },
+        { 0x1.9c2d163a1aa2dp-1, 0x1.bc2860d22477p-3 },
+        { 0x1.886e6037841edp-1, 0x1.1058bc8a07ee1p-2 },
+        { 0x1.767dcf5534862p-1, 0x1.4043057b6ee09p-2 }
+    };
+
+    double z, r, r2, y, y0, invc, logc;
+    UINT32 ix, iz, tmp;
+    int k, i;
+
+    ix = *(UINT32*)&x;
+    /* Fix sign of zero with downward rounding when x==1. */
+    if (ix == 0x3f800000)
+        return 0;
+    if (ix - 0x00800000 >= 0x7f800000 - 0x00800000) {
+        /* x < 0x1p-126 or inf or nan. */
+        if (ix * 2 == 0)
+            return math_error(_SING, "logf", x, 0, -1.0 / x);
+        if (ix == 0x7f800000) /* log(inf) == inf. */
+            return x;
+        if (ix * 2 > 0xff000000)
+            return x;
+        if (ix & 0x80000000)
+            return math_error(_DOMAIN, "logf", x, 0, (x - x) / (x - x));
+        /* x is subnormal, normalize it. */
+        x *= 0x1p23f;
+        ix = *(UINT32*)&x;
+        ix -= 23 << 23;
+    }
+
+    /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
+       The range is split into N subintervals.
+       The ith subinterval contains z and c is near its center. */
+    tmp = ix - 0x3f330000;
+    i = (tmp >> (23 - 4)) % (1 << 4);
+    k = (INT32)tmp >> 23; /* arithmetic shift */
+    iz = ix - (tmp & 0xff << 23);
+    invc = T[i].invc;
+    logc = T[i].logc;
+    z = *(float*)&iz;
+
+    /* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */
+    r = z * invc - 1;
+    y0 = logc + (double)k * Ln2;
+
+    /* Pipelined polynomial evaluation to approximate log1p(r). */
+    r2 = r * r;
+    y = A[1] * r + A[2];
+    y = A[0] * r2 + y;
+    y = y * r2 + (y0 + r);
+    return y;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

@@ -120,14 +120,6 @@ static float CDECL unix_lgammaf(float x)
 #endif
 }
 
-/*********************************************************************
- *      logf
- */
-static float CDECL unix_logf( float x )
-{
-    return logf( x );
-}
-
 /*********************************************************************
  *      log10
  */
@@ -243,7 +235,6 @@ static const struct unix_funcs funcs =
     unix_fmaf,
     unix_lgamma,
     unix_lgammaf,
-    unix_logf,
     unix_log10,
     unix_log10f,
     unix_log1p,

dlls/msvcrt/unixlib.h

@@ -30,7 +30,6 @@ struct unix_funcs
     float           (CDECL *fmaf)(float x, float y, float z);
     double          (CDECL *lgamma)(double x);
     float           (CDECL *lgammaf)(float x);
-    float           (CDECL *logf)(float x);
     double          (CDECL *log10)(double x);
     float           (CDECL *log10f)(float x);
     double          (CDECL *log1p)(double x);