msvcrt: Import log10f implementation from musl.

dlls/msvcrt/math.c

@@ -1295,10 +1295,62 @@ float CDECL logf( float x )
  */
 float CDECL log10f( float x )
 {
-  float ret = unix_funcs->log10f( x );
-  if (x < 0.0) return math_error(_DOMAIN, "log10f", x, 0, ret);
-  if (x == 0.0) return math_error(_SING, "log10f", x, 0, ret);
-  return ret;
+    static const float ivln10hi = 4.3432617188e-01,
+        ivln10lo = -3.1689971365e-05,
+        log10_2hi = 3.0102920532e-01,
+        log10_2lo = 7.9034151668e-07,
+        Lg1 = 0xaaaaaa.0p-24,
+        Lg2 = 0xccce13.0p-25,
+        Lg3 = 0x91e9ee.0p-25,
+        Lg4 = 0xf89e26.0p-26;
+
+    union {float f; UINT32 i;} u = {x};
+    float hfsq, f, s, z, R, w, t1, t2, dk, hi, lo;
+    UINT32 ix;
+    int k;
+
+    ix = u.i;
+    k = 0;
+    if (ix < 0x00800000 || ix >> 31) { /* x < 2**-126 */
+        if (ix << 1 == 0)
+            return math_error(_SING, "log10f", x, 0, -1 / (x * x));
+        if ((ix & ~(1u << 31)) > 0x7f800000)
+            return x;
+        if (ix >> 31)
+            return math_error(_DOMAIN, "log10f", x, 0, (x - x) / (x - x));
+        /* subnormal number, scale up x */
+        k -= 25;
+        x *= 0x1p25f;
+        u.f = x;
+        ix = u.i;
+    } else if (ix >= 0x7f800000) {
+        return x;
+    } else if (ix == 0x3f800000)
+        return 0;
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    ix += 0x3f800000 - 0x3f3504f3;
+    k += (int)(ix >> 23) - 0x7f;
+    ix = (ix & 0x007fffff) + 0x3f3504f3;
+    u.i = ix;
+    x = u.f;
+
+    f = x - 1.0f;
+    s = f / (2.0f + f);
+    z = s * s;
+    w = z * z;
+    t1= w * (Lg2 + w * Lg4);
+    t2= z * (Lg1 + w * Lg3);
+    R = t2 + t1;
+    hfsq = 0.5f * f * f;
+
+    hi = f - hfsq;
+    u.f = hi;
+    u.i &= 0xfffff000;
+    hi = u.f;
+    lo = f - hi - hfsq + s * (hfsq + R);
+    dk = k;
+    return dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi + hi * ivln10hi + dk * log10_2hi;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

@@ -121,14 +121,6 @@ static float CDECL unix_lgammaf(float x)
 }
 
 /*********************************************************************
- *      log10f
- */
-static float CDECL unix_log10f( float x )
-{
-    return log10f( x );
-}
-
-/*********************************************************************
  *      log2
  */
 static double CDECL unix_log2(double x)
@@ -203,7 +195,6 @@ static const struct unix_funcs funcs =
     unix_fmaf,
     unix_lgamma,
     unix_lgammaf,
-    unix_log10f,
     unix_log2,
     unix_log2f,
     unix_pow,

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 *log10f)(float x);
     double          (CDECL *log2)(double x);
     float           (CDECL *log2f)(float x);
     double          (CDECL *pow)(double x, double y);