msvcrt: Import log10 implementation from musl.

dlls/msvcrt/math.c

@@ -2964,10 +2964,82 @@ double CDECL log( double x )
  */
 double CDECL log10( double x )
 {
-  double ret = unix_funcs->log10( x );
-  if (x < 0.0) return math_error(_DOMAIN, "log10", x, 0, ret);
-  if (x == 0.0) return math_error(_SING, "log10", x, 0, ret);
-  return ret;
+    static const double ivln10hi = 4.34294481878168880939e-01,
+        ivln10lo = 2.50829467116452752298e-11,
+        log10_2hi = 3.01029995663611771306e-01,
+        log10_2lo = 3.69423907715893078616e-13,
+        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, s, z, R, w, t1, t2, dk, y, hi, lo, val_hi, val_lo;
+    UINT32 hx;
+    int k;
+
+    hx = u.i >> 32;
+    k = 0;
+    if (hx < 0x00100000 || hx >> 31) {
+        if (u.i << 1 == 0)
+            return math_error(_SING, "log10", x, 0, -1 / (x * x));
+        if ((u.i & ~(1ULL << 63)) > 0x7ff0000000000000ULL)
+            return x;
+        if (hx >> 31)
+            return math_error(_DOMAIN, "log10", x, 0, (x - x) / (x - x));
+        /* subnormal number, scale x up */
+        k -= 54;
+        x *= 0x1p54;
+        u.f = x;
+        hx = u.i >> 32;
+    } else if (hx >= 0x7ff00000) {
+        return x;
+    } else if (hx == 0x3ff00000 && u.i<<32 == 0)
+        return 0;
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    hx += 0x3ff00000 - 0x3fe6a09e;
+    k += (int)(hx >> 20) - 0x3ff;
+    hx = (hx & 0x000fffff) + 0x3fe6a09e;
+    u.i = (UINT64)hx << 32 | (u.i & 0xffffffff);
+    x = u.f;
+
+    f = x - 1.0;
+    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;
+
+    /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
+    hi = f - hfsq;
+    u.f = hi;
+    u.i &= (UINT64)-1 << 32;
+    hi = u.f;
+    lo = f - hi - hfsq + s * (hfsq + R);
+
+    /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */
+    val_hi = hi * ivln10hi;
+    dk = k;
+    y = dk * log10_2hi;
+    val_lo = dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi;
+
+    /*
+     * Extra precision in for adding y is not strictly needed
+     * since there is no very large cancellation near x = sqrt(2) or
+     * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
+     * with some parallelism and it reduces the error for many args.
+     */
+    w = y + val_hi;
+    val_lo += (y - w) + val_hi;
+    val_hi = w;
+
+    return val_lo + val_hi;
 }
 
 /*********************************************************************

dlls/msvcrt/unixlib.c

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

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